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Abstract 

The  relationship  between  depth-averaged  velocity  and  bottom  stress  for  wind-driven  flow  in 
unstratified  coastal  waters  is  examined  here.  The  adequacy  of  traditional  linear  and  quadratic 
drag  laws  is  addressed  by  comparison  with  a  2|-D  model.  A  2|-D  model  is  one  in  which  a 
simplified  1-D  depth-resolving  model  (DRM)  is  used  to  provide  an  estimate  of  the  relationship 
between  the  flow  and  bottom  stress  at  each  grid  point  of  a  depth-averaged  model  (DAM). 
Bottom  stress  information  is  passed  from  the  DRM  to  the  DAM  in  the  form  of  drag  tensor 
with  two  components:  one  which  scales  the  flow  and  one  which  rotates  it.  This  eliminates 
the  problem  of  traditional  drag  laws  requiring  the  flow  and  bottom  stress  to  be  collinear.  In 
addition,  the  drag  tensor  field  is  updated  periodically  so  that  the  relationship  between  the 
velocity  and  bottom  stress  can  be  time-dependent.  However,  simplifications  in  the  2|-D  model 
that  render  it  computationally  efficient  also  impose  restrictions  on  the  time-scale  of  resolvable 
processes.  Basically,  they  must  be  much  longer  than  the  vertical  diffusion  time  scale. 

Four  progressively  more  complicated  scenarios  are  investigated.  The  important  factors 
governing  the  importance  of  bottom  friction  in  each  are  found  to  be  1)  non-dimensional  surface 
Ekman  depth,  u*s/fh  where  ums  is  the  surface  shear  velocity,  /  is  the  Coriolis  parameter  and 
h  is  the  water  depth  2)  the  non-dimensional  bottom  roughness,  zo/h  where  z0  is  the  roughness 
length  and  3)  the  angle  between  the  wind  stress  and  the  shoreline.  Each  has  significant  influence 
on  the  drag  law.  The  drag  tensor  magnitude,  r,  and  the  drag  tensor  angle,  0  are  functions  of 
all  three,  while  a  drag  tensor  which  scales  with  the  square  of  the  depth-averaged  velocity  has  a 
magnitude,  Cd,  that  only  depends  on  z0/h. 

The  choice  of  drag  law  is  found  to  significantly  affect  the  response  of  a  domain.  Spin 
up  times  and  phase  relationships  vary  between  models.  In  general,  the  2|-D  model  responds 
more  quickly  than  either  a  constant  r  or  constant  Cd  model.  Steady-state  responses  are  also 
affected.  The  two  most  significant  results  are  that  failure  to  account  for  0  in  the  drag  law 
sometimes  leads  to  substantial  errors  in  estimating  the  sea  surface  height  and  to  extremely 
poor  resolution  of  cross-shore  bottom  stress.  The  latter  implies  that  cross-shore  near-bottom 
transport  is  essentially  neglected  by  traditional  DAMs. 
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Chapter  1 


Introduction 


This  thesis  presents  a  method  for  improving  wind-driven  depth-averaged  circulation  models 
(DAMs)  by  providing  a  more  realistic  estimate  of  bottom  stress  than  traditional  drag  for¬ 
mulations.  Improved  DAMs  should  yield  better  predictions  of  coastal  water  levels  and  mass 
transports.  In  its  existing  form,  the  modelling  approach  used  here  is  applicable  to  unstratified, 
wind-driven  coastal  flows.  However,  it  can,  with  minimal  effort,  be  adopted  also  to  include 
tidally-driven  flows. 

Chapter  1  contains  a  brief  motivating  discussion  based  on  the  present  state  of  wind-driven 
circulation  modelling.  In  addition,  a  basic  description  of  the  modelling  procedure  and  an  outline 
of  the  thesis  structure  axe  given. 

1.1  Motivation 

The  two  most  relevant  problems  to  depth-averaged  circulation  modelling  are  storm  surge  and 
horizontal  transport  prediction.  Each  has  obvious  health  and  economic  implications.  Therefore, 
any  method  by  which  they  can  be  improved  seems  worthwhile  pursuing. 

Storm  surges  are  among  the  most  devastating  of  natural  disasters.  Thousands  of  people  have 
been  killed  and  billions  of  dollars  of  property  destroyed  by  storm  surges  in  the  United  States 
alone  (Murty,  1984).  Similarly,  Japan,  India,  Bangladesh  and  all  of  the  European  countries 
surrounding  the  North  Sea  have  experienced  catastrophic  storm  surges.  In  fact,  few  of  the 
world’s  coastlines  have  been  spared. 
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Consequently,  a  great  deal  of  interest  in  providing  useful  storm  surge  forecasting  schemes 
has  long  existed  in  the  oceanographic  community.  This  is  illustrated  clearly  by  the  fact  that 
in  1957  enough  literature  on  the  subject  had  accumulated  to  prompt  Charnock  and  Crease 
to  write  a  review  paper  on  storm  surge  observations  and  prediction.  Since  then,  storm  surge 
research  has  flourished  (See  additional  reviews  by  Welander,  1961;  Groen  and  Groves,  1962; 
Bretschneider,  1967;  Heaps  1967,  1983;  and  a  recent  900  page  monograph  by  Murty,  1984). 

Perhaps  the  biggest  trend  in  storm  surge  forecasting  over  the  last  30  years  has  been  the  shift 
from  empirical  to  numerical  prediction  schemes.  The  move  gained  a  great  deal  of  momentum 
in  the  late  1960s  with  the  introduction  of  a  number  of  wind-driven  coastal  DAMs  (e.g.  Heaps, 
1969).  Three-dimensional  (3-D)  models  soon  followed  and  gained  limited  popularity  in  the  mid- 
1970s  (e.g.  Heaps  and  Jones,  1975).  Both  types  have  since  developed  significantly.  Nonetheless, 
DAMs  constitute  the  vast  majority  of  operational  storm  surge  forecasting  models  and  are  not 
likely  to  be  replaced  by  3-D  models  in  the  near  future.  This  is  due  primarily  to  the  latter’s 
extreme  computational  expense  (Davies,  1989  personal  communication). 

Wind-driven  transport  modelling  does  not  have  nearly  the  long  history  that  storm  surge 
modelling  does,  nor  has  it  enjoyed  as  much  success.  There  are  many  reasons  for  this.  Perhaps 
chief  among  them  is  the  difficulty  of  making  Lagrangian  transport  calculations  due  to  their 
extreme  sensitivity  to  small  scale  structure  in  the  velocity  field.  Sea  surface  heights,  on  the 
other  hand,  are  a  more  robust  variable  and,  therefore,  more  easily  predicted. 

It  has  only  been  recently  that  widespread  availability  of  supercomputers  has  led  to  improve¬ 
ments  in  small  scale  depth-averaged  flow  structure  resolution  (e.g.  Signell  and  Geyer,  1989). 
In  addition,  the  introduction  of  efficient  coordinate  systems  (e.g.  Blumberg  and  Herring,  1983) 
and  other  numerical  tools  has  aided  this  effort. 

Superimposed  on  the  horizontal  resolution  problem  is  the  difficulty  that  DAMs  are  often 
used  to  predict  the  transport  of  quantities  that  are  not  vertically  weU-mixed.  This  is  usually 
attempted  by  introducing  an  artificial  horizontal  dispersion  coefficient  which  is  to  account  for 
vertical  differences  in  transport  (e.g.  Kossik  et  a!.,  1987). 

Despite  their  inherent  problems,  depth-averaged  transport  models  stand  to  be  improved  in 
two  ways  by  better  estimates  of  bottom  stress.  The  first  is  simply  through  improvement  of 
the  predicted  flow  structure.  If  the  velocity  field  is  modelled  more  accurately,  the  associated 
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transport  fields  should  be  modelled  more  accurately  too.  The  second  is  relevant  to  models 
where  the  transported  quantity  is  known  to  be  concentrated  near  the  sea  bottom.  In  such 
cases,  more  accurate  predictions  of  bottom  stress  can  be  translated  into  better  predictions  of 
the  advecting  near-bottom  flow  field. 

There  is  actually  a  third  way  that  the  method  described  here  might  be  useful  for  transport 
modelling.  However,  it  is  not  specifically  addressed  in  this  thesis.  The  method  yields  simplified 
estimates  of  the  vertical  profile  of  horizontal  velocity  at  each  grid  point  in  a  DAM.  This  informa¬ 
tion  could  be  used  to  calculate  better  horizontal  dispersion  coefficients  or  to  advect  quantities 
differentially  at  different  depths  in  an  associated  transport  model  or  to  obtain  estimates  of 
momentum  transport  coefficients  for  the  DAM  itself. 

Clearly,  the  above  discussion  is  not  an  exhaustive  review  of  storm  surge  and  depth-averaged 
transport  modelling.  Instead,  it  is  intended  to  present  the  context  in  which  improvements 
of  bottom  stress  estimates  may  be  useful.  The  comments  above  are  not  entirely  speculation, 
however.  In  fact,  Blumberg  and  Oey  (1985)  in  a  recent  review  cite  lack  of  appropriate  bottom 
stress  formulations  for  DAMs  as  one  of  the  most  significant  deficiencies  in  coastal  circulation 
modelling  at  present. 

1.2  Problem  Definition 

When  the  momentum  and  mass  conservation  equations  governing  wind-driven  flow  dynamics 
are  integrated  vertically,  equations  describing  evolution  of  the  sea  surface  elevation  and  depth- 
averaged  velocity  result.  These  equations  contain  terms  representing  both  momentum  sources 
and  sinks.  The  sources  are  wind  stress  and  the  atmospheric  pressure  gradient.  The  sink  is  drag 
on  the  flow  imposed  by  the  seabed.  Specification  of  the  meteorological  source  terms  is,  in  itself, 
a  very  difficult  problem  and  is  fax  from  being  satisfactorily  resolved  (Murty,  1984;  Blumberg 
and  Oey,  1985).  Nonetheless,  this  thesis  focuses  only  on  the  problem  of  specifying  the  bottom 
stress. 

Because  the  bottom  stress  is  a  function  of  the  flow,  it  must  somehow  be  related  to  the 
dependent  variables,  depth-averaged  velocity  and  sea  surface  height,  and  to  the  known  meteo¬ 
rological  forcing.  This  presents  a  severe  problem  as  the  bottom  stress  is  a  function  of  the  near 
bottom  flow  field.  Unfortunately,  the  bottom  flow  field  is  related  to  the  dependent  variables  and 
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Figure  1-1:  Schematic  illustrating  potential  problems  with  “traditional”  bottom  stress  formu¬ 
lations 

forcing  through  a  very  complicated  set  of  physics.  The  crux  of  the  bottom  stress  formulation 
problem  is  the  parameterization  of  the  relationship  between  the  bottom  stress  and  the  near 
bottom  flow,  dependent  variable  and  forcing  fields. 


1.3  “Traditional”  Drag  Laws 

The  most  common  method  of  parameterizing  bottom  stress  in  DAMs  is  the  quadratic  drag  law. 
In  such  a  law,  the  bottom  stress  is  equal  to  the  product  of  the  water  density,  the  depth-averaged 
velocity  vector,  its  magnitude,  and  a  prescribed  nondimensional  drag  coefficient.  Quadratic 
drag  laws  are  easy  to  implement  and  add  little  computational  effort  to  the  solution  procedure. 
However,  ease  of  use  is  potentially  overshadowed  by  an  oversimplification  of  physics. 

A  good  example  of  problems  caused  by  the  physical  oversimplification  of  quadratic  drag 
laws  is  illustrated  in  figure  1-la.  The  depth-averaged  flow  is  equal  to  0  but  there  is  a  finite  flow 
near  the  bottom.  The  quadratic  drag  law  predicts  no  bottom  stress  in  this  case  even  though 
it  exists.  A  more  important  example  of  oversimplification  is  illustrated  in  figure  1-lb.  Coriolis 
acceleration  causes  the  flow  to  rotate  with  depth.  Thus,  the  direction  of  the  bottom  stress  is 
not  that  of  the  depth-averaged  flow.  The  quadratic  drag  law  predicts  the  wrong  bottom  stress 
by  restricting  it  to  be  in  the  same  direction  as  the  depth-averaged  velocity. 
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Another  traditional  formulation  used  in  place  of  the  quadratic  drag  law  is  the  linear  drag 
law.  In  a  linear  drag  law,  bottom  stress  is  formulated  as  the  product  of  water  density,  depth- 
averaged  velocity  and  an  a  priori  specified  resistance  coefficient  with  dimensions  [length/time]. 
The  resistance  coefficient  can  be  thought  of  as  a  drag  coefficient  multiplied  by  a  representative 
depth-averaged  current  speed  since  it  essentially  replaces  the  product  of  the  drag  coefficient  and 
the  variable  depth-averaged  velocity  magnitude  in  the  quadratic  law.  Linear  drag  laws  retain 
the  problems  illustrated  in  figure  1-1. 

Because  the  resistance  coefficient  is  not  temporally-  or  spatially-varying  like  the  depth- 
averaged  velocity  magnitude,  linear  formulations  are  often  considered  more  restrictive  than 
quadratic  drag  laws.  Consequently,  they  tend  to  be  used  more  often  in  process-oriented  models 
(e.g.  Chapman,  1985)  while  quadratic  laws  are  more  frequently  used  in  simulation  models  (e.g. 
Beardsley  and  Haidvogel,  1981). 

This  is  not  to  say  that  linear  drag  laws  are  used  only  in  process  models.  In  fact,  certain 
characteristics  of  linear  drag  laws  make  them  “preferable”  to  quadratic  drag  laws.  If  the  depth- 
averaged  governing  equations  to  be  modelled  are  otherwise  linear,  using  a  linear  drag  law  keeps 
them  that  way.  This  often  improves  the  numerical  stability  of  the  associated  model. 

Unfortunately,  linear  drag  laws  are  also  chosen  sometimes  when  a  modeller  simply  assumes 
that  bottom  stress  cannot  be  predicted  accurately  by  either  of  the  traditional  drag  laws.  The 
computationally  simpler  linear  drag  law  is  chosen  and  the  resistance  coefficient  is  used  as  a 
tunable  parameter  to  improve  agreement  between  model  results  and  observations. 

1.4  Alternatives  to  Traditional  Drag  Laws 

Various  methods  have  been  implemented  to  overcome  the  problems  of  traditional  drag  laws. 
Fixes  range  from  adding  a  constant  vector  to  the  bottom  stress  to  resolving  the  near  bottom  flow 
field  with  a  3-D  model.  Clearly,  there  is  a  trade-off  between  computational  ease  and  realistic 
representation  of  bottom  stress. 

The  difficulty  with  traditional  drag  laws  illustrated  in  figure  1-la  was  addressed  by  Bowden 
(1953)  and  Reid  (1957).  Both  examined  the  situation  of  steady  flow  forced  by  a  wind  perpen¬ 
dicular  to  a  vertical  wall.  In  the  depth-averaged  sense,  there  is  no  flow  perpendicular  to  the 
wall.  Therefor^,  traditional  drag  laws  predict  no  bottom  stress  in  that  direction.  However,  in 
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reality,  a  strong  wind  blowing  directly  perpendicular  to  a  wall  would  force  a  surface  flow  toward 
the  wall  which  would  be  balanced  by  a  return  flow  near  the  bottom. 

Bowden  used  a  constant  eddy  viscosity  model  to  resolve  the  flow  field  and  to  show  that  the 
bottom  stress  (force  on  the  bottom  by  the  flow)  for  laminar  flow  should  be  between  0%  and  50% 
of  the  surface  stress  and  directed  upwind.  Reid  used  a  more  complicated  mixing  length  model 
to  suggest  that  for  turbulent  flow  the  bottom  stress  should  be  around  10%  of  the  surface  stress 
and  directed  upwind.  These  ideas  were  generalized  by  others  (e.g.  Groen  and  Groves,  1962) 
who  suggested  that  traditional  drag  laws  should  be  modified  by  adding  a  constant  fraction 
(typically  10%,  after  Reid)  of  the  wind  stress  to  the  bottom  stress.  This  idea  breaks  down 
immediately,  however,  for  a  strong  wind  stress  parallel  to  a  wall  (or  a  straight  coastline).  In 
such  a  case,  the  flow  is  entirely  downwind  throughout  the  water  column  and  the  stress  exerted 
by  the  flow  on  the  seabed  is  downwind  rather  than  upwind. 

Jenter  and  Madsen  (1989)  used  a  steady-state  bi-linear  eddy  viscosity  model  (a  modified 
version  of  which  is  described  in  Chapter  2)  to  examine  the  idea  of  adding  a  constant  fraction 
of  the  wind  stress  to  the  bottom  stress.  They  found  that  the  ratio  of  cross-shore  bottom  stress 
to  cross-shore  wind  stress  was  a  complicated  function  of  bottom  roughness,  wind  strength  and 
angle  between  the  wind  and  the  wall.  Using  a  different  eddy  viscosity,  Davies  (1987)  also  found 
a  large  amount  of  variability  in  the  fractional  relationship  between  cross-shore  bottom  stress 
and  surface  stress.  Jamart  and  Ozer  (1987)  analytically  examined  steady-state  wind-driven  flow 
in  a  closed  basin  and  showed  that  the  relationship  between  bottom  stress  and  surface  stress 
necessarily  included  other  terms  beside  the  traditional  drag  term  and  a  fraction  of  the  wind 
stress.  There  is  ample  evidence  that  using  a  predetermined  constant  fraction  of  the  wind  stress 
to  improve  bottom  stress  in  DAMs  is  inappropriate. 

A  slightly  more  involved  approach  than  the  drag  laws  described  above  is  used  in  so-called 
2  ^-dimensional  (2  |-D)  models  where  bottom  stress  is  predicted  for  a  DAM  by  estimating  the 
bottom  flow  structure  with  a  separate  one-dimensional  (1-D)  model.  The  upper  complexity 
limit  for  2  |-D  models  might  be  considered  a  fully  3-D  model.  However,  as  used  in  this  thesis, 
2  |-D  model  refers  to  a  DAM  coupled  to  a  simplified  depth-resolving  model  (DRM).  Therefore, 
the  term  implies  that  the  combined  model  improves  the  bottom  stress  estimates  by  resolving 
the  near  bottom  flow  structure  while  not  entirely  sacrificing  the  computational  efficiency  of  a 
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DAM. 


Many  2  |-D  models  have  their  origin  in  the  work  of  Welander  (1957),  who  showed  that  for  a 
constant  eddy  viscosity,  the  bottom  stress  can  be  related  to  time  histories  of  the  wind  stress  and 
sea  surface  slope  through  a  convolution  integral.  Since  the  wind  stress  is  known  and  the  surface 
slope  at  any  point  can  be  calculated  from  the  depth-averaged  equations,  the  bottom  stress 
integral  can  be  evaluated  in  parallel  with  the  DAM  yielding  improved  bottom  stress  estimates. 
Jelesnianski  (1970)  provided  a  method  of  analytically  simplifying  the  convolution  integral  and 
applied  it  in  a  storm  surge  model.  Forristall  (1974)  developed  a  similar  convolution  integral 
and  applied  it  in  a  study  of  storm-generated  currents.  Kielmann  and  Kowalik  (1980)  used  a 
combination  of  finite  difference  and  analytic  methods  to  approximate  the  convolution  integral. 

The  convolution  integral  models  necessarily  employ  a  constant  eddy  viscosity.  A  more  flexi¬ 
ble  structure  is  obviously  desirable.  To  this  end,  Nihoul  (1977)  introduced  a  technique  by  which 
a  temporally  and  vertically-variable  eddy  viscosity  can  be  used,  provided  the  eddy  viscosity  is 
separable  in  depth  and  time.  Separability  allowed  Nihoul  to  specify  the  velocity  profile  as 
the  sum  of  a  set  of  vertical  modes  with  time-dependent  amplitudes  but  time-invariant  vertical 
structure.  The  shape  functions  for  the  modes  are  chosen  to  be  eigenfunctions  of  a  second-order 
ordinary  differential  equation  involving  the  depth-dependent  part  of  the  eddy  viscosity.  Much 
like  evaluating  the  convolution  integrals,  calculating  the  mode  amplitudes  corresponds  to  com¬ 
puting  an  inverse  Laplace  transform  involving  the  sea  surface  slope  and  boundary  shear  stresses. 
Resultant  expressions  for  the  near  bottom  flow  can  be  inverted  numerically  with  some  difficulty 
to  relate  the  bottom  stress  to  the  depth-averaged  velocity,  sea  surface  slope  and  wind-stress. 

Independently,  Jordan  and  Baker  (1980)  and  Baker  and  Jordan  (1980,  1981)  developed 
similar  expansion  techniques.  Davies  (1987)  applied  the  method  and  showed  that,  with  the 
appropriate  problem  formulation,  mode  amplitudes  could  be  calculated  by  recursion  relation¬ 
ships  greatly  simplifying  determination  of  bottom  stress  for  the  DAM.  Davies  also  claimed  that 
only  2  to  4  modes  were  necessary  to  yield  DAM  results  comparable  to  those  of  a  fully  3-D 
model.  Unfortunately,  his  test  case  was  in  water  so  deep  that  bottom  stress  was  probably  not 
dynamically  important.  Because  of  this,  major  differences  in  bottom  stress  formulations  might 
not  produce  noticeable  differences  in  model  output,  rendering  his  results  inconclusive. 

The  convolution  and  eigenfunction  methods  are  a  clear  improvement  over  traditional  drag 
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laws,  but  have  two  potential  problems.  The  first,  and  most  obvious,  is  that  the  form  of  eddy 
viscosity  is  restricted  to  be,  at  best,  separable  in  depth  and  time  and,  at  worst,  constant.  Com¬ 
plicated  flow  regimes  in  shallow  wind-driven  waters  necessitate  violation  of  these  restrictions. 
The  presence  of  both  surface  and  bottom  boundary  layers,  which  interact  and  vary  in  relative 
strength,  requires  the  eddy  viscosity  to  adjust  its  relative  strength  near  the  surface  and  bottom 
boundaries  as  a  function  of  time,  which  cannot  happen  in  the  convolution  and  eigenfunction 
formulations. 

The  second  problem  is  more  subtle.  The  convolution  and  eigenfunction  methods  require 
that  either  a  slip  or  no-slip  boundary  condition  be  applied  at  the  sea  bed.  A  slip  condition 
relates  near  bottom  velocity  to  near  bottom  shear,  and  a  no  slip  condition  requires  the  near 
bottom  velocity  to  vanish.  Each  can  impose  unrealistic  solutions  on  the  near  bottom  flow  profile 
when  applied  individually.  More  realistically,  the  bottom  flow  should  satisfy  both  no-slip  and 
shear  conditions.  In  order  to  do  this,  however,  the  eddy  viscosity  must  be  a  function  of  the 
flow.  Such  a  relationship  precludes  a  priori  knowledge  of  the  eddy  viscosity  function. 

Unfortunately,  the  problem  of  allowing  eddy  viscosity  to  be  flow-dependent  is  perhaps  best 
addressed  by  second-order  closure  models  (e.g.  Mellor  and  Yamada,  1974)  in  which  time- 
and  depth-dependent  differential  equations  for  turbulent  velocity  and  length  scales  axe  solved. 
Models  have  been  formulated  with  second-order  closure  schemes  (e.g.  Blumberg  and  Mellor, 
1987  and  Johns  and  Oguz,  1987),  but  at  a  substantial  computational  cost.  Noyes  and  Stevens 
(1987),  in  fact,  claim  that  the  cost  of  including  length  and  velocity  scale  equations  increases  3-D 
model  computational  and  storage  requirements  by  a  minimum  of  30%.  Scholtz  et  a!.  (1987) 
have  formulated  a  model  with  second-order  closure  that  might  be  construed  as  an  extreme  2 
1-D  model  because  the  differential  equations  for  the  length  and  velocity  scales  are  replaced 
by  algebraic  approximations.  Nonetheless,  the  required  computational  effort  must  be  at  least 
comparable  to  a  3-D  model  without  second-order  closure  since  the  momentum  equations  are 
not  simplified. 


1.5  The  Present  Approach 

The  goal  of  the  work  reported  in  this  thesis  is  to  produce  a  2  |-D  formulation  that  incorporates 
realistic  yet  simple  turbulent  closure  assumptions.  The  DRM  employs  an  eddy  viscosity  closure 
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that  yields  quasi- analytical  velocity  profiles  satisfying  the  2  ^-D  model  criterion  for  computa¬ 
tional  efficiency.  Despite  permitting  an  analytical  solution,  the  eddy  viscosity  is  quite  realistic. 
It  has  a  time-dependent  functional  form  which  strongly  relates  it  to  the  near  bottom  flow  field. 
Consequently,  it  can  satisfy  both  shear  and  no-slip  bottom  boundary  conditions.  Also,  because 
the  problems  of  interest  are  wind-driven,  the  eddy  viscosity  accounts  for  the  presence  of  a  sur¬ 
face  boundary  layer  and  for  the  possibility  that  it  will  interact  with  the  bottom  boundary  layer 
in  shallow  enough  water. 

The  DRM  communicates  bottom  stress  information  to  the  DAM  by  estimating  a  drag  tensor 
at  each  DAM  grid  point.  The  tensor  is  used  in  the  DAM  to  scale  and  rotate  the  depth-averaged 
velocity  to  yield  an  improved  estimate  of  bottom  stress.  This  eliminates  the  problem  illustrated 
in  figure  1-lb  because  the  depth-averaged  velocity  and  bottom  stress  do  not  have  to  be  in  the 
same  direction.  However,  the  problem  described  in  1-la  is  not  eliminated  since  there  is  still  no 
bottom  stress  when  the  depth-averaged  velocity  is  exactly  equal  to  0. 

The  drag  tensor  formulation  is  partially  chosen  for  reasons  of  computational  efficiency.  With 
it,  the  bottom  stress  need  not  be  updated  at  every  time  step  of  the  DAM.  Bottom  stress  can 
still  vary  between  updates  because  it  is  also  a  function  of  the  depth-averaged  velocity. 

1.6  Organization  of  the  Thesis 

Modelling  details  are  described  in  Chapter  2  including  discussion  of  the  DAM,  the  DRM  and 
the  combined  2  |-D  model.  Chapter  3  contains  2  |-D  model  results  for  a  variety  of  wind-driven 
coastal  flow  situations.  Chapter  4  is  a  summary  of  these  results  and  their  implications  for 
wind-driven  depth-averaged  flow  modelling. 
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Chapter  2 


Model  Formulation 


2.1  The  Depth-averaged  Model 


Since  the  goal  of  this  thesis  is  to  provide  an  improved  method  of  estimating  bottom  stress  in 
wind-driven  DAMs,  a  standard  DAM  was  chosen  as  a  starting  point.  The  model  was  developed 
by  Beardsley  and  Haidvogel  (1981)  to  model  the  storm  response  of  the  Mid- Atlantic  Bight.  It 
has  since  been  employed  by  Chapman  (1985)  to  test  a  variety  of  open  boundary  conditions. 
Well  documented  model  application  was  a  primary  factor  in  its  selection  here. 

The  governing  equations  are  the  vertically-integrated,  linearized  momentum  and  continuity 
equations  for  an  unstratified  fluid.  They  are 
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where  (  is  sea  surface  displacement,  U  alongshelf  (x)  transport  ([17]  =  length2/time),  V  cross¬ 
shelf  (y)  transport,  g  gravitational  acceleration,  h  water  depth,  /  the  Coriolis  parameter,  p  water 
density,  and  (r^,  r^)  and  [rax,Tay)  are  the  (x,y)  bottom  and  surface  stresses,  respectively. 

Numerical  considerations  notwithstanding,  and  presuming  the  wind  stress  is  known,  solution 
of  the  system  (2.1),  (2.2)  and  (2.3)  depends  entirely  on  correct  specification  of  the  a  priori 
unknown  bottom  stress.  As  discussed  in  Chapter  1,  this  is  usually  done  by  explicitly  relating 
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the  bottom  stress  to  the  transport  through  a  drag  law.  In  a  recent  review  of  numerical  coastal 
circulation  modeling,  Blumberg  and  Oey  (1985)  consider  the  specification  of  appropriate  drag 
laws  to  be  one  of  the  single  biggest  problems  facing  depth-averaged  modellers  at  present. 

The  two  most  common  bottom  stress  formulations  in  DAMs  are  drag  coefficients  and  re¬ 
sistance  coefficients.  In  a  drag  coefficient  formulation,  bottom  stress  is  related  to  transport 
through  the  quadratic  expression 


n  =  pCd 


u\u 


h 2 


(2.4). 


where  the  drag  coefficient,  Cd,  is  dimensionless.  In  a  resistance  coefficient  formulation,  bottom 
stress  is  related  to  transport  through  the  linear  expression 


(2.5) 


where  the  resistance  coefficient,  r,  has  dimensions  [length/time].  Both  (2.4)  and  (2.5)  imply 
the  very  restrictive  assumption  that  bottom  stress  and  transport  are  in  the  same  direction. 
However,  due  to  the  possible  importance  of  Coriolis  accelerations  in  the  momentum  balance, 
the  two  are  not  likely  to  be  colinear  (i.e.  ffe  X  U  0)  because  bottom  stress  is  determined  by 
the  near-bottom  flow  structure,  and  not  an  average  over  the  water  column. 

Adopting  the  basic  idea  that  bottom  stress  and  transport  need  not  be  colinear,  the  first 
improvement  suggested  in  this  thesis  is  the  introduction  of  a  drag  tensor.  With  a  drag  tensor, 
both  scaling  and  rotation  of  the  transport  vector  are  performed  in  order  to  produce  a  bottom 
stress  vector.  For  the  drag  coefficient,  this  is  written  as 
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To  emphasize  the  scaling  and  rotational  nature  of  (2.6)  and  (2.7),  they  can  also  be  written  as 
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Figure  2-1:  Definition  sketch  for  the  angular  difference,  0,  between  bottom  stress,  and 
transport,  U.  Also  the  angle  between  the  bottom  stress  and  x-axis,  a,  and  between  the  transport 
and  x-axis,  /?  are  shown. 
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Taking  the  vector  norm  of  (2.8)  and  (2.9)  and  rearranging  the  results  yields 
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and 
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which  imply  that  Cd  and  r  play  the  same  role  in  (2.8)  and  (2.9)  as  they  do  in  traditional  drag 
laws.  In  other  words,  they  are  still  the  factors  by  which  the  magnitude  of  the  depth-averaged 
velocity  is  scaled  to  yield  the  magnitude  of  the  bottom  stress. 

The  drag  tensor  angle,  0,  also  has  physical  significance  and  an  analog  in  traditional  drag 
laws.  It  is  the  angle  between  the  transport  and  bottom  stress  (measured  counter-clockwise  from 
the  transport,  as  shown  in  figure  2-1).  Traditional  drag  laws  simply  assume  0  =  0.  Notice  that 
choosing  0  =  0  reduces  (2.8)  and  (2.9)  to  (2.4)  and  (2.5). 

The  amount  of  improvement  offered  by  a  drag  tensor  formulation  lies  in  the  extent  to  which 
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Cd  or  r  and  0  can  be  correctly  specified.  A  primary  goal  of  this  thesis  is  to  find  a  computationally 
simple  way  of  accurately  specifying  drag  tensors  in  typical  wind-driven  coastal  environments. 
Introducing  (2.9)  into  (2.1)  and  (2.2)  yields  modified  depth-integrated  momentum  equations1 


w  ~fV  =  ^hf,  +  Tf-ThUcos>+riiVsine 
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It  is  interesting  to  note  that  (2.12)  and  (2.13)  can  also  be  written  as 
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In  other  words,  the  transport  equations  with  a  drag  tensor  can  be  arranged  in  such  a  way  that 
they  appear  identical  to  the  transport  equations  with  a  colinear  resistance  law.  This  means  that 
the  drag  tensor  (2.9)  merely  appears  as  a  redefinition  of  the  Coriolis  parameter  and  resistance 
coefficient.  Consequently,  the  numerical  scheme  employed  by  Chapman  (1985)  need  not  be 
changed  at  all  to  accomodate  a  drag  tensor,  so  long  as  the  new  definitions  are  accepted. 

The  finite  difference  equations  corresponding  to  (2.14),  (2.15)  and  (2.3)  are 
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1For  brevity,  parallel  derivations  of  equations  using  (2.8)  are  no  longer  carried  out.  Instead  the  linear  formu¬ 
lation  is  exclusively  employed.  Similar  derivations  for  the  quadratic  formulation  are  straightforward.  However, 
since  r  is  a  time- varying  quantity  in  the  2j-D  model,  the  choice  of  a  linear  instead  of  quadratic  drag  tensor  is 
insignificant. 
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Figure  2-2:  The  Arakawa  “C”  grid 

The  superscript  n  implies  that  the  variable  is  evaluated  at  time,  t  =  nAt  where  At  is  the  model 
time  step.  The  subscripts  i  and  j  imply  that  the  variable  is  evaluated  at  the  element  in  the 
x-direction  and  the  element  in  the  y-direction  of  the  grid.  The  grid  spacing  in  both  the  x 
and  y  directions  is  assumed  constant  and  equal  to  Ax.  An  overbar  indicates  that  an  average 
of  the  variable  over  adjacent  grid  elements  is  required  (see  Beardsley  and  Haidvogel  1981  for 
details.). 

The  grid  used  for  the  DAM  is  an  Arakawa  “C”  grid  shown  in  figure  2-2.  The  “C”  grid 
is  a  spatially  staggered  grid  chosen  for  its  gravity  wave  propagation  characteristics,  its  ability 
to  handle  coastal  boundary  conditions  simply  and  the  convenience  with  which  centered  spa¬ 
tial  differences  are  calculated.  Grid  characteristics  are  thoroughly  discussed  by  Mesinger  and 
Arakawa  (1976). 

One  difficulty  encountered  with  the  staggered  grid  is  its  inability  to  handle  the  Coriolis  terms 
simply  in  the  transport  equations  because  U  and  V  are  not  defined  at  the  same  nodes.  Thus, 
a  spatial  averaging  of  transports  for  the  Coriolis  terms  is  required  as  indicated  in  (2.18)  and 
(2.19).  The  DAM  uses  an  averaging  scheme  developed  by  Platzman  (1972)  where  surrounding 
transports  are  first  weighted  by  inverse  depth  then  averaged.  The  complete  expression  is  given 
in  the  original  model  description  by  Beardsley  and  Haidvogel  (1981). 

Time  differencing  is  done  explicitly  by  a  forward/backward  scheme.  In  such  a  scheme, 
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transports  are  stepped  forward  alternately  with  surface  elevation  as  if  staggered  in  time  by  At/2. 
Von  Neumann  stability  analysis  of  the  forward/backward  scheme  for  gravity  wave  propagation 
yields  a  Courant-Friedrichs-Lewy  (CFL)  stability  criterion  (Mesinger  and  Arakawa,  1976) 


At  < 


1  Ax 

a/2  y/gh 


(2.21) 


Note  that  the  on-axis  friction  terms,  fU /h,  are  formulated  implicitly  with  respect  to  the  trans¬ 
port.  Partial  justification  for  doing  so  comes  from  Von  Neumann  stability  analysis  of  (2.18), 
(2.19)  and  (2.20)  with  /  =  0,  =  0  and  f  =  constant,  which  shows  that  an  explicit  form 

yields  a  stability  criterion  At  <  2 h/f  in  addition  to  (2.21)  (Chapman,  1987  personal  communi¬ 
cation).  With  implicit  friction  terms,  however,  the  numerical  scheme  is  unconditionally  stable 
with  respect  to  constant  r.  The  implicit  scheme  is  adopted  here  with  the  hope  that  the  stability 
characteristics  for  constant  f  are  inherited  by  the  2|-D  model  with  variable  f. 

So  far,  nothing  has  been  said  about  actual  values  of  the  drag  tensor.  Jenter  and  Madsen 
(1989)  used  a  DRM  identical  to  the  one  to  be  described  in  section  2.2  in  order  to  compute  values 
of  r  and  0  for  the  simple  case  of  steady  wind-driven  flow  in  an  ocean  of  infinite  horizontal  extent 
(i.e.  no  pressure  gradients).  Figures  2-3  and  2-4  depict  drag  tensor  variation  as  a  function  of 
dimensionless  surface  Ekman  depth,  u*3/fh ,  where  ums  is  the  surface  shear  velocity  defined 
u*s  =  \/W¥p-  The  different  curves  in  each  figure  are  for  different  ratios  of  bottom  roughness 
length  to  water  depth,  zo/h. 

In  order  to  facilitate  interpretation  of  figures  2-3  and  2-4,  Table  2.1  contains  values  of 
u+s/fh  for  different  wind  speeds  and  water  depths.  The  wind  speeds  are  converted  into  stresses 
by  application  of  Wu’s  (1982)  formula  and  /  is  taken  to  be  10-4  s-1,  as  it  is  throughout  this 
thesis.  The  cases  in  the  upper  righthand  half  of  the  table  (i.e.  above  the  thick  line)  are  of 
little  interest  or  dynamical  importance  for  storm  surge  and  transport  modelling.  Therefore, 
this  thesis  concentrates  on  cases  taken  from  the  remainder  of  Table  2.1. 

Table  2.2  contains  typical  values  of  zo/h  for  different  sediment  types.  The  roughness  lengths 
are  based  on  Nikuradse  sand  grain  roughnesses,  zq  =  d/30  where  d  is  the  grain  diameter. 
It  may  seem  strange  that  pebbles  and  cobbles  are  included  in  the  table,  since  they  are  not 
typical  bottom  types  for  coastal  areas.  Grant  and  Madsen  (1979,  1982),  showed,  however,  that 
the  presence  of  surface  gravity  waves  or  bed  forms  can  cause  the  flow  to  “feel”  an  enhanced 
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Figure  2-4:  Variation  of  drag  tensor  direction  as  a  function  of  Ekman  depth  to  water  depth 
ratio  (from  Jenter  and  Madsen,  1989). 
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Wind  speed 
at  10m, 
m/s 

Beaufort 

number, 

description 

Water  depth,  m 

2 

5 

10 

20 

50 

100 

200 

1 

1 

Light  air 

5 

2 

1 

0.5 

0.2 

0.1 

0.1 

3 

2 

Light  breeze 

16 

7 

X 

2 

0.7 

0.3 

0.2 

5 

3 

Gentle  breeze 

29 

12 

V  3 

1 

0.6 

0.3 

8 

4 

Moderate  breeze 

51 

20 

10 

2 

1 

0.5 

!  11 

5 

Fresh  breeze 

75 

30 

15 

\ 

.  3 

1 

0.8 

14 

6 

Strong  breeze 

101 

41 

20 

10 

V 

2 

1 

17 

7 

Moderate  gale 

130 

52 

26 

13 

3 

1 

21 

8 

Fresh  gale 

171 

68 

34 

17 

\ 

3 

2 

24 

9 

Strong  gale 

204 

82 

41 

20 

s\ 

V  4 

2 

28 

10 

Whole  gale 

250 

100 

50 

25 

10 

V 

3 

33 

11 

Storm 

313 

125 

63 

31 

13 

Table  2.1:  Dimensionless  surface  Ekman  depth,  u„j  fh,  as  a  function  of  wind  speed  10m  above 
the  water  surface,  and  water  depth,  h.  Wind  speeds  are  converted  to  stresses  using  Wu’s  (1982) 
formula  and  /  =  10~4  s-1. 

roughness  comparable  to  cobbles.  Therefore,  the  larger  values  of  zo/h  are  still  relevant  to  the 
wind-driven  coastal  flows  considered  here. 

Figure  2-3  shows  that  the  steady-state  value  of  r  can  vary  over  multiple  orders  of  magni¬ 
tude  as  a  function  of  water  depth  and  wind  stress.  This  implies  that  a  priori  knowledge  of 
r  is  highly  dependent  on  a  priori  knowledge  of  u»».  It  also  suggests  that  situations  involving 
spatially  variable  winds  or  topography  should  require  spatially  variable  values  of  r  in  order  to 


Grain 

diameter, 

cm 

Size  class 

Water  depth, 

m 

2 

20 

200 

Cobbles 

ECTTiBl 

Pebbles 

EH- 

0.1 

Coarse  Sand 

2  X  10~5 

2  X  10~6 

2  x  lO"7 

0.05 

Medium  Sand 

8  X  10~6 

8  x  10"7 

8  X  10“8 

Fine  Sand 

4  x  10“6 

4  x  10~7 

4  x  10“8 

KiSIll 

Silt 

2  X  lO"7 

2  x  10“8 

2  X  10~9 

Table  2.2:  Range  of  dimensionless  bottom  roughness  length,  zo/h,  as  a  function  of  grain  size, 
d ,  and  water  depth,  h. 
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model  bottom  stress  correctly.  Similarly,  situations  with  temporally  variable  winds  will  require 
temporally  variable  values  of  r. 

Like  figure  2-3,  figure  2-4  implies  a  great  deal  of  drag  tensor  variation  with  wind  speed  and 
water  depth.  It  is  apparent  that,  for  a  steady-state  infinite  horizontal  extent  ocean,  a  colinear 
drag  law  (i.e.  0  =  0)  applies  only  in  extremely  shallow  water  or  for  extremely  high  winds.  As 
with  r,  it  appears  that  a  priori  knowledge  of  0  is  highly  dependent  on  a  priori  knowledge  of 
u*s.  Consequently,  any  spatial  or  temporal  variation  in  the  wind  or  water  depth  should  require 
variation  in  0. 

At  least  for  an  ocean  of  infinite  horizontal  extent,  spatial  and  temporal  variation  of  6  is 
of  more  concern  than  variation  in  r  because  much  of  the  variability  in  r  can  be  removed  by 
assuming  a  quadratic  drag  law.  Figure  2-5  (also  from  Jenter  and  Madsen,  1989)  shows  variation 
of  Cd  as  a  function  of  Ekman  depth  to  water  depth  ratio.  Clearly,  Cd  is  not  nearly  as  sensitive 
as  r  to  u*a/fh.  In  fact,  variation  in  Cd  over  the  important  range  of  u„a/ fh  >  5  is  comparable 
to  that  associated  with  a  factor  of  2  uncertainty  in  the  bottom  roughness. 

Having  established  that  the  drag  tensor  is  sensitive  to  forcing  and  topographic  variation, 
a  method  by  which  r  and  0  can  be  systematically  varied  in  space  and  time  is  introduced.  A 
DRM  is  used  to  update  periodically  the  drag  tensor  at  each  grid  element  of  the  DAM.  The 
combination  of  the  DAM  and  DRM  constitutes  the  2|-D  model. 

Notice  that  adoption  of  the  2|-D  modelling  strategy  means  that  the  choice  to  use  r  instead 
of  Cd  is  now  immaterial  because  the  drag  tensor  is  allowed  to  vary  in  time  and  space.  Therefore, 
provided  it  is  updated  often  enough,  temporal  changes  in  forcing  should  produce  appropriate 
changes  in  the  drag  tensor. 

2.2  The  Depth-Resolving  Model 

This  section  describes  a  quasi-analytical  depth- resolving  model,  DRM,  used  for  predicting  verti¬ 
cal  profiles  of  horizontal  velocity  in  a  wind-  and  pressure-driven  unstratified  fluid.  The  resulting 
profiles  are,  in  turn,  used  to  determine  a  drag  tensor  for  the  DAM. 

Derivation  of  the  DRM  is  based  on  a  set  of  linearized  momentum  equations  that  can  be 
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(2.22) 

(2.23) 


vertically  integrated  to  obtain  (2.1)  and  (2.2).  Specifically, 

du  .  3(  3/  du\ 

m-fv  =  -Sdi  +  Tz[^a;) 

d  (  dv\ 

Tz  VTTz) 


dv  ,  ,  d(  ^ 

m+fu  =  ~9di  + 


where  (u(z),v(z))  is  the  depth-dependent  horizontal  velocity  in  the  (x,y)  direction.  Notice  that 
the  stresses  are  modelled  as  the  product  of  an  eddy  viscosity,  vt,  and  the  vertical  shear  in  the 
velocity  profile.  Also,  the  hydrostatic  assumption  is  made  to  relate  the  pressure  gradient  to  the 
sea  surface  slope. 

The  key  to  resolving  successfully  the  flow  structure  using  (2.22)  and  (2.23)  lies  in  choosing 
an  appropriate  form  for  the  eddy  viscosity.  Since  the  2^-D  model  is  intended  to  describe  wind- 
driven  flow  in  shallow  water,  the  eddy  viscosity  must  account  for  the  presence  of  both  surface 
and  bottom  boundary  layers  as  well  as  for  the  possibility  that  the  two  overlap  in  sufficiently 
shallow  water. 

The  DRM  assumes 


vt  =  \ 


KU,s(h  -  z), 
KU,bZ, 


Zm  <  Z  <  h 
0  <  z  <  zm 


(2.24) 


where  k  =  0.4  is  von  Karman’s  constant,  u„&  is  the  a  priori  unknown  bottom  shear  velocity 
defined  by  |f&| /p  =  ulb  and  zm  is  a  matching  height  (figure  2-6).  The  origin  of  the  2- axis  is 
chosen  to  coincide  with  the  sea  bottom  as  shown  in  figure  2-6.  The  matching  height  is  defined 
by 


Zm 


«.*  +  u„6 


(2.25) 


a  formulation  which  was  first  suggested  by  Madsen  (1977). 

It  is  important  to  note  that  the  matching  height  is  not  the  point  where  the  two  functions 
of  (2.24)  overlap.  Such  a  choice  would  apportion  a  larger  part  of  the  water  column  to  the 
function  having  the  smaller  shear  velocity.  The  choice  of  matching  height  (2.25)  ensures  that 
each  function  is  apportioned  a  part  of  the  water  column  according  to  the  relative  strength  of  its 
shear  velocity.  Consequently,  the  eddy  viscosity  is  discontinuous  at  z  =  zm  (figure  2-6)  unless 
=  u**.  Further  implications  of  this  discontinuity  are  discussed  below. 

The  bilinear  character  of  (2.24)  is  chosen  for  two  reasons.  First,  in  keeping  with  the  effi¬ 
ciency  goal  of  the  2|-D  model,  it  does  not  preclude  a  simple  analytical  solution  to  (2.22)  and 
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(2.23).  As  important  as  simplicity,  however,  is  its  ability  to  reproduce  much  of  the  surface 
and  bottom  boundary  layer  behavior  observed  in  real  flows.  For  example,  wind-driven  surface 
boundary  layers  in  the  laboratory  have  been  shown  to  exhibit  a  log-deficit  profile  near  the  free 
surface  (Shemdin,  1972;  Wu,  1975).  Log-deficit  layers  have  also  been  observed  in  the  field  us¬ 
ing  near  surface  drifters  (Churchill  and  Csanady,  1983;  Csanady,  1984),  current  meter  drifters 
(Richman  et  al.,  1987)  and  platform  mounted  acoustic  current  meters  (Terray  et  al.,  1987).  An 
approximation  to  log-deficit  variation  is  produced  by  an  eddy  viscosity  which  increases  linearly 
away  from  the  free  surface  while  the  stress  remains  relatively  constant. 

Another  characteristic  of  surface  boundary  layers  which  is  dependent  upon  the  eddy  viscos¬ 
ity  is  the  rotational  behavior.  Field  experiments  in  deep  water  indicate  a  structure  which  is 
somewhat  different  from  the  classical  Ekman  spiral.  One  major  difference  is  the  angle  which 
the  surface  velocity  makes  with  the  wind  stress.  Classical  theory  (Ekman,  1905)  predicts  that 
this  angle  should  be  45°,  while  observations  (e.g.  Stacey  et  al.,  1986)  indicate  that  it  may  be 
substantially  smaller.  An  eddy  viscosity  increasing  linearly  away  from  the  free  surface  has  been 
shown  by  Madsen  (1977)  to  predict  a  reduction  in  angle  between  the  wind  and  surface  velocity. 

Bottom  boundary  layer  behavior  has  also  been  studied  in  both  laboratory  (e.g.  Nezu  and 
Rodi,  1986)  and  field  experiments  (see  Soulsby,  1983;  Grant  and  Madsen,  1986  for  reviews.). 
For  steady  unstratified  flow,  the  turbulent  boundary  layer  near  a  solid  wall  is  well  known  to 
contain  a  region  of  logarithmic  variation  with  depth  (Clauser,  1956).  As  with  the  surface 
log-deficit  layer,  this  variation  is  approximated  fairly  accurately  by  an  eddy  viscosity  which 
increases  linearly  away  from  the  boundary  while  the  stress  profile  remains  relatively  constant. 
In  addition,  the  slope  of  the  eddy  viscosity  has  been  shown  to  be  proportional  to  the  bottom 
shear  velocity,  u„f,. 

Constant  stress  layers  are  clearly  an  approximation  near  the  boundaries.  Nonetheless,  (2.24) 
reproduces  the  basic  near-boundary  flow  features.  There  is  evidence  (e.g.  Businger  and  Arya, 
1974)  that  the  linear  eddy  viscosity  should  be  cut  off  at  a  distance,  7 «*//,  away  from  the 
boundary  where  typical  7  values  have  been  quoted  in  the  range  0.04  to  0.06  (Lavelle  and 
Mofjeld,  1983).  The  fact  that  the  present  analysis  neglects  the  eddy  viscosity  cutoff  implies 
that  the  resulting  theory  should  be  limited  to  shallow  water  only. 

A  working  definition  of  “shallow  water”  is  suggested  by  the  results  of  Jenter  and  Madsen 
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(1989)  who  found  that  the  cutoff  points  of  the  surface  and  bottom  boundary  layers  overlapped 
(i.e.  7 («»s  +  u»j,)//  >  h )  for  «*,/ fh  >  17  when  7  =  0.04  and  for  u*s/fh  >  11  when  7  =  0.06. 
Further  experimentation  with  a  trilinear  eddy  viscosity  (linearly  connecting  the  cutoff  points 
when  they  do  not  overlap)  showed  that  the  bilinear  solutions  were  quantitatively  robust  above 
u*s/fh  «  5  for  7  =  0.06.  Recent  results  obtained  by  Mofjeld  (1988)  for  purely  pressure- 
driven  flow  in  deep  water  (i.e.  u*s  =  0  in  (2.24))  indicate  that  differences  in  bottom  stress 
magnitude  with  or  without  a  cutoff  are  at  most  10%  for  all  reasonable  flow  conditions.  Bottom 
stress  directional  variation  is  similarly  small.  Consequently,  concern  over  application  of  the 
2|-D  model  presented  here  should  only  be  raised  for  “light  winds  or  deep  water”  below  say 
U.s/fh  ss  5  (see  table  2.1.),  which  is  outside  the  range  of  primary  interest. 

The  question  of  the  applicability  of  (2.24)  in  the  bottom  stress  prediction  process  becomes 
moot  for  low  values  of  u*,/ fh  because  bottom  stress  does  not  play  an  important  dynamical 
role  in  the  wind-driven  DAM.  Csanady  (1984)  uses  simple  scaling  arguments  to  support  the 
idea  that  bottom  stress  loses  its  dynamical  importance  when  u*4//h  =  0(1).  Likewise,  the 
fairly  elaborate  turbulence  closure  model  of  Overland  et  ai.  (1984)  predicts  that  wind-driven 
ice  movement  shows  no  bottom  friction  effects  when  «*«/ fh  is  less  than  about  10. 

Having  defined  the  eddy  viscosity  and  its  range  of  applicability,  a  means  of  solving  (2.22)  and 
(2.23)  is  required.  One  potential  choice  is  by  Laplace  transform  methods.  However,  this  is  only 
possible  for  eddy  viscosity  profiles  of  particular  forms.  Jelesnianski  (1970)  and  Forristall  (1974) 
presented  Laplace  transform-like  solutions  relating  bottom  stress  to  the  time  histories  of  local 
wind  stress  and  sea  surface  slope  for  wind-driven  coastal  flow  with  constant  vj.  Madsen  (1977) 
presented  a  Laplace  transform  solution  of  (2.22)  and  (2.23)  in  an  infinitely  deep  ocean  with 
a  time-independent  eddy  viscosity  that  increased  linearly  away  from  the  sea  surface  (i.e.  the 
surface  function  of  (2.24)).  Unfortunately,  these  Laplace  transform  solutions  are  not  applicable 
in  the  present  study  because  the  functional  dependence  of  v t  on  the  a  priori  unknown  bottom 
stress  implies  that  vt  is  unpredictably  time-dependent. 

An  alternative  is  the  decomposition  of  the  velocity  profile  into  vertical  modes  that  are 
separable  in  depth  and  time.  This  procedure  was  first  introduced  by  Heaps  (1972)  for  a  constant 
eddy  viscosity.  It  was  later  modified  for  time-  and  depth-dependent  eddy  viscosity  by  both 
Nihoul  (1977),  Baker  and  Jordan  (1980,  1981)  and  Jordan  and  Baker  (1980).  The  procedure 
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has  since  been  used  by  a  variety  of  investigators  (e.g.  Davies,  1987).  Introduction  of  a  modal 
structure  reduces  the  full  solution  of  (2.22)  and  (2.23)  to  the  solution  of  first-order  ordinary 
differential  equations  for  the  time-dependent  mode  amplitudes.  However,  this  is  only  possible 
when  the  eddy  viscosity  is  formulated  as  a  separable  function  of  depth  and  time.  As  with 
the  Laplace  transform  solution,  functional  dependence  of  vr  on  the  time- varying  bottom  stress 
makes  the  modal  approach  inapplicable  because  v-j-  is  not  separable  when  zm  is  allowed  to  vary 
(i.e.  the  ratio  of  bottom  stress  to  surface  stress  is  allowed  to  vary). 

Since  Laplace  transform  and  modal  decomposition  solutions  are  not  possible,  (2.22)  and 
(2.23)  must  be  simplified  if  an  analytic  solution  is  sought.  One  way  to  do  this  is  by  assuming 
that  the  local  time  derivatives  in  (2.22)  and  (2.23)  are  depth-independent  and  equal  to  their 
depth-averaged  values  (i.e.  replace  du/dt  in  (2.22)  and  (2.23)  by  (l/h)dU /dt).  This  permits 
(2.22)  and  (2.23)  to  be  rewritten  as 


where 


,,  ^  d  f  d(u-ug)\ 

-/(»-»„)  =  — ) 

/v  x  9  (  d(V~V9)\ 

/(«-%)  =  ai  — si — ) 

i/i  dv  d(\ 

u°  = 

1  fldu  d(\ 

V°  ~  f  \h  dt  +  3  dx) 


(2.26) 

(2.27) 

(2.28) 
(2.29) 


The  subscript  g  refers  to  the  fact  that  ug  and  va  without  the  local  time  derivatives  are  merely 
geostrophic  velocities.  Equations  (2.28)  and  (2.29)  emphasize  that  the  unsteady  terms  are 
treated  as  a  modification  to  the  pressure  gradient  with  their  effect  being  felt  uniformly  through¬ 
out  the  water  column.  Limitations  on  temporally  resolvable  processes  caused  by  this  approxi¬ 
mation  are  discussed  later. 

Equations  (2.26)  and  (2.27)  are  rewritten  by  taking  their  dot  product  with  (1,  *)  (i  =  V^T) 
to  yield 

Tz  ("T^r) _  ifyv  =  0  (2-30) 

where  W  =  (u  —  ug)  +  i(v  —  vg)  is  the  complex  difference  between  the  total  velocity  and  (ug,  vg). 
Inserting  (2.24)  into  (2.30)  yields 


Tz  (KU”(h  ~  2)7r)  ” ifw  =  °- 


<  z  <  h 


(2.31) 
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^  -  *'/W  =  0,  0  <  z  <  z, 


(2.32) 


Solutions  to  (2.31)  and  (2.32)  are  found  in  many  applied  mathematics  texts  (e.g.  Hildebrand, 
1976).  They  are 


W=  < 


A+  (ber 

2  \/f(h  -  z)/ku,s 

+  ibei 

2 \Jf{h  -  z)/ku*s 

B+  (ker 

2v//(/i  -  z)/Kuma 

4-  ikei 

2 y/f(h  -  z)/nu*a 

A  (ber 

2s/fzfKunh  +  ibei  2yl/fz/numb\)  + 

B  (ker[  2%/JzJku^  +ikei  [2i//z/ku*&]  ) 


zm  <  z  <  h 

0  <  z  <  zm 


(2.33) 


where  A±  and  B±  are  complex  constants  determined  by  surface  and  bottom  boundary  condi¬ 
tions  and  ber,bei,ker,  and  kei  are  zero-order  Kelvin  functions. 

The  boundary  conditions  for  determining  A*  and  B±  are 


W\z=*o  =  “(%  +  iva) 

(2.34) 

..  d\V  tsx  +  irsy 

lim  vt—t—  = - 

z->h  dz  p 

(2.35) 

where  the  complex  forcing  quantities,  ( ug  +  ivg)  and  (rsx  +  trsy),  are  passed  to  the  DEM  from 
the  DAM. 

Equation  (2.34)  is  merely  a  no-slip  condition  applied  a  small  distance,  zq,  above  the  sea 
bottom.  This  is  necessary  because  (2.33)  behaves  logarithmically  as  z  — ►  0.  In  reality,  z0,  the 
roughness  length,  is  a  function  of  sediment  characteristics,  bed  configuration,  and  near-bottom 
turbulence  structure  (e.g.,  due  to  wave-current  interaction).  However,  the  analysis  presented 
here  assumes  that  zq  is  known  a  priori,  thereby  ignoring  details  of  near-bottom  flow  structure 
and  bottom  roughness  in  the  determination  of  roughness  length. 

Condition  (2.35)  is  a  statement  relating  known  surface  shear  stress  to  unknown  surface 
shear.  It  is  necessary  to  take  the  limit  as  the  boundary  is  approached  since  the  velocity  profile 
behaves  logarithmically  and  thus  the  shear  profile  varies  as  l/(/i  —  z ).  The  stress  limit  is  finite, 
however,  since  vt  approaches  0  in  such  a  way  that  it  balances  the  shear  as  it  approaches  oo. 

Supplementing  (2.34)  and  (2.35)  are  matching  conditions  at  the  level  of  discontinuity  in  the 
eddy  viscosity. 

W|z=2m+  =  W\z=Zm_  (2.36) 
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Uj 


dW 


dz 


—  Vt 


dW 


z=zm+ 


dz 


(2.37) 


z — zm— 


Matching  conditions  (2.36)  and  (2.37),  respectively,  ensure  continuity  of  velocity  and  shear 
stress  across  zm.  Because  the  shear  stress  is  continuous  and  the  eddy  viscosity  is  discontinuous 
at  zm  the  shear  must  be  discontinuous  there.  This  means  the  velocity  profile  is  kinked  at  zm. 
Clearly,  a  kink  is  unrealistic.  However,  it  is  not  pronounced  so  long  a szm  >  zq,  which  is 
normally  the  case.  Comparison  between  bilinear  eddy  viscosity  results  and  those  for  a  con¬ 
tinuous  trilinear  eddy  viscosity  mentioned  earlier  show  little  effect  on  the  velocity  profile,  its 


depth-average  or  the  associated  bottom  stress. 

While  boundary  conditions  (2.34),  (2.35)  and  matching  conditions  (2.36),  and  (2.37)  de¬ 
termine  A±  and  B±,  the  problem  is  still  not  closed2.  The  bottom  shear  velocity,  u+b,  is  a 
priori  unknown  and  part  of  the  solution.  Therefore,  the  solution  must  be  iterated  to  satisfy  an 
additional  compatibility  condition 


lim  vt 

z->  0 


dW 

dz 


I  Tbx  +  *Tby  I  _  „,2 
~  —  u*b 


(2.38) 


Satisfaction  of  (2.38)  ensures  that  the  bottom  stress  magnitude  on  which  the  eddy  viscosity 
is  based  is  the  actual  bottom  stress  magnitude  produced  by  the  current  profile.  Unlike  most 
vertical  profile  models,  both  no-slip  (2.34)  and  shear  (2.38)  boundary  conditions  are  satisfied. 

The  idea  of  using  a  turbulent  mixing  scheme  based  on  an  a  priori  unknown  bottom  shear 
velocity  is  not  new.  Thomas  (1975)  employed  a  similar  assumption,  but  did  not  incorporate  a 
surface  mixing  formulation.  Thomas’  analysis  stressed  the  need  for  an  improved  eddy  viscosity 
formulation  that  accounts  for  surface  boundary  layer  as  well  as  overlapping  boundary  layer 
effects  in  shallow  wind-driven  waters.  In  addition  to  Thomas’  analysis,  similar  iterative  solution 
procedures  have  been  employed  for  boundary  layer  models  of  wave-current  interaction  (e.g. 
Smith  and  Long,  1976;  Smith,  1977;  and  Grant  and  Madsen,  1979). 

Figure  2-7  shows  a  schematic  representation  of  the  solution  procedure  for  the  DRM  as  it  is 
employed  in  the  2|-D  model.  The  iterative  nature  of  the  solution  is  apparent.  This  reduces 
the  computational  efficiency  of  the  DRM,  but  experience  has  shown  only  5  or  6  iterations  are 
required  to  obtain  agreement  in  to  less  than  0.1%.  Upon  completion  of  iteration,  the  bottom 


2  Complete  expressions  for 


A±  and  B*  are  given  by  Jenter  (1987). 
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stress  and  velocity  profile  (and,  consequently,  the  transport3)  are  obtained.  The  bottom  stress 
and  transport  are  then  manipulated  to  yield  estimates  of  r  and  0  to  be  used  in  the  DAM. 

Before  describing  the  2|-D  model,  a  brief  discussion  of  the  obvious  numerical  alternative  to 
an  analytic  DRM  should  be  presented.  Clearly,  a  numerical  procedure  is  adaptable  to  the  eddy 
viscosity  (2.24),  and  eliminates  the  problems  with  Laplace  transforms  and  modal  decomposition. 
It  can  also  handle  depth-dependence  of  the  local  time  derivatives  giving  it  an  advantage  over  the 
Kelvin  function  solution.  However,  use  of  a  numerical  procedure  in  the  context  of  a  2^-D  model 
raises  the  question  of  its  efficiency.  To  examine  this  issue,  steady-state  (i.e.  dropping  du/dt ) 
solutions  of  (2.22)  and  (2.23)  were  calculated  analytically  and  numerically  and  the  results 
compared  for  speed  and  accuracy.  The  numerical  solution  consisted  of  a  centered  difference 
scheme  on  a  logarithmically  stretched  grid  with  fine  spacing  near  the  surface  and  bottom  in 
order  to  resolve  the  boundary  layers  without  restrictively  small  grid  spacing  in  the  interior  of 
the  water  column.  The  problem  was  formulated  as  a  boundary  value  problem  and  solved  by  a 
standard  tri-diagonal  matrix  inversion.  Depth-averaging  of  the  resultant  solution  was  done  by 
Simpson’s  rule.  Additional  details  of  the  numerical  procedure  are  given  by  Jenter  (1987). 

Table  2.3  summarizes  the  comparison.  One  hundred  runs  with  random  wind  stresses,  pres¬ 
sure  gradients,  bottom  roughnesses  and  water  depths  were  performed  and  statistics  of  the 
results  calculated.  The  table  shows  that  for  the  numerical  and  analytical  solutions  to  have  sim¬ 
ilar  computation  times,  the  numerical  solution  can  utilize  only  25  grid  points.  When  this  is  the 
case,  however,  the  relative  accuracy  of  the  numerical  scheme  is  terrible.  While  the  numerical 
and  analytical  bottom  stress  magnitudes  agree  in  the  mean,  two  standard  deviations  represent 
±53%.  This  can  be  reduced  to  ±8%  at  a  cost  of  4  times  the  CPU  time  and  to  ±2%  at  20  times 
the  CPU  time.  Given  these  results  and  the  fact  that  the  DRM  is  the  most  CPU-intensive  part 
of  the  2|-D  model,  a  numerical  solution  to  (2.22)  and  (2.23)  is  rejected  on  efficiency  grounds. 

2.3  The  2|-D  Model 

Together  the  DAM  and  DRM  provide  a  wind-driven  coastal  circulation  model  with  an  improved 
estimate  of  bottom  stress.  The  basic  2|-D  modelling  scheme  is  illustrated  by  the  following 

3  A  complete  analytic  expression  for  the  transport  is  given  by  Jenter  (1987). 
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case 

CPU  time 

IAI 

^analytical 

M _ 

^analytical 

Analytical 

0:02.6 

1 

1 

Numerical  25pts 

0:02.4 

1.04  ±  .53 

1.00  ±  .18 

(.77-2.7) 

(.85-1.45) 

Numerical  lOOpts 

0:08.9 

1.00  ±  .08 

1.00  ±  .04 

(.94-1.30) 

(.97-1.10) 

Numerical  500pts 

0:43.2 

1.00  ±  .02 

1.00  db  .008 

(.99-1.06) 

(.99-1.02) 

Table  2.3:  Comparison  of  speed  and  accuracy  for  analytical  and  numerical  solutions  with  ran¬ 
dom  forcings.  Statistics  are  mean ± two  standard  deviations  and  (range)  for  100  runs. 

simple  procedure: 

•  step  the  DAM  through  time  using  a  time  step  based  on  the  CFL  stability  condition  until 
a  decision  is  made  to  update  the  drag  tensor 

•  pass  forcing  information  from  the  DAM  to  the  DRM  so  that  drag  tensor  estimates  can 
be  obtained  at  the  grid  points  throughout  the  domain 

•  pass  the  drag  tensor  estimates  from  the  DRM  back  to  the  DAM  and  return  to  the  first 
step  above 

This  process  is  illustrated  in  more  detail  in  figure  2-8. 

The  2|-D  model  actually  consists  of  three  parts.  The  DAM  and  DRM  have  been  described 
in  sections  2.1  and  2.2,  respectively,  while  the  third  part,  the  roughness  model,  is  mentioned 
merely  for  completeness.  In  reality,  it  is  only  a  formalism  at  present  since  the  roughness  field  is 
specified  a  priori  and  left  unchanged.  Nonetheless,  the  model  is  formulated  in  such  a  way  that 
the  roughness  field  can  be  updated  whenever  the  DRM  is  run.  This  is  done  with  an  eye  toward 
the  possibility  of  including  such  roughness  altering  processes  as  wave-current  interaction  (e.g. 
Grant  and  Madsen,  1979)  and  moveable  bed  effects  (e.g.  Grant  and  Madsen,  1982)  in  future 
versions. 

The  first  step  in  running  the  2--D  model  is  initialization  of  the  field  variables,  U,  V,  (,  r  and 
0.  All  of  these  are  set  equal  to  0  when  the  model  is  spun  up  from  rest.  In  addition,  the  physical 
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Figure  2-8:  Schematic  of  the  full  2|-D  model. 
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Figure  2-9:  The  Arakawa  “C”  grid  showing  superimposed  r,6  nodes. 


configuration  of  the  model  must  be  input,  which  includes  an  Arakawa  “C”  grid  superimposed 
over  the  desired  topography,  a  bottom  roughness  field  and  an  initial  wind  stress  field. 

After  initialization,  the  DAM  is  time  stepped  forward  using  the  forward/backward  scheme. 
Following  a  predetermined  number  of  time  steps,  the  DRM  is  called  to  update  the  r  and  0 
fields.  The  update  is  performed  at  the  unoccupied  fourth  corner  of  each  grid  element  (figure  2- 
9)  because  calculation  of  the  pressure  gradient  needed  to  force  the  DRM  is  straightforward  using 
the  surrounding  (  points  as  is  calculation  of  the  local  time  derivatives  of  transport.  Specifically, 
at  the  r,6  node  of  the  grid  element  i,j, 

~  Ct-lJ  +  Cij—l  ~  Ct-IJ-l;  ,  P,j  ~  Ctj-i  +  C«-i,j  ~  C«— l,j— 1  ; 

2A*  +  2Ax  J  (2^9) 


and 


dU  U"j  -  u?r 
dt  2mA< 


+  +  -  v™ + vr-u  -  yr~ 


— m 

i  J  • 


2mAt 


(2.40) 


where  (i,  j)  are  unit  vectors  in  the  (i,  y)  directions  and  m  is  the  number  of  time  steps  between 
drag  tensor  updates.  Clearly  (2.39)  and  (2.40)  cannot  be  applied  at  the  DAM  boundaries. 
Therefore,  V£,  dU/dt ,  r  and  9  are  calculated  only  at  interior  grid  elements. 

Once  fields  of  V£  and  dU/dt  (thus  ug  and  vg  from  (2.28)  and  (2.29))  are  calculated,  the 
feedback  from  the  DAM  to  the  DRM  at  each  node  is  specified.  The  DRM  is  forced  by  {ug  +  ivg) 
along  with  the  same  wind  stress  and  topography  as  used  in  the  DAM.  After  importing  the 
forcing,  the  roughness  field,  which  is  kept  constant  here,  is  input. 
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Given  the  forcing,  topography  and  roughness  at  each  interior  r,  0  node,  the  DRM  is  applied 
subject  to  the  boundary  conditions  (2.34),  (2.35),  (2.36),  (2.37),  and  (2.38).  After  convergence, 
the  bottom  stress  and  velocity  profile  are  calculated  from  (2.38)  without  the  absolute  value  signs 
and  (2.33),  respectively.  The  velocity  profile  is  then  integrated  to  yield  a  transport  vector. 

The  bottom  stress  and  transport  predicted  by  the  DRM  are  transformed  into  estimates  of 
r  and  0  using  definitions  derived  from  (2.11)  and  figure  2-1.  Since  forcings  are  only  supplied 
at  interior  r,0  nodes,  the  DRM  results  are  extrapolated  to  border  nodes.  Extrapolation  to  the 
boundary  nodes  is  done  linearly  from  the  two  nearest  r,  0  points  along  a  line  intersecting  the 
node  and  perpendicular  to  the  boundary.  Corner  extrapolations  take  into  account  the  presence 
of  both  boundaries.  After  extrapolation,  the  complete  field  of  drag  tensor  values  is  passed  back 
to  the  DAM.  Time  stepping  of  the  DAM  continues  until  it  is  necessary  to  update  the  drag 
tensor  field  again. 

It  should  be  noted  that  the  DAM  requires  r  and  0  information  at  U  and  V  nodes  while 
the  DRM  provides  it  at  r,  0  nodes.  Therefore,  an  intermediate  spatial  averaging  takes  place  in 
order  to  provide  drag  tensor  values  at  the  proper  locations.  At  the  Uij  node: 


(/  -  f)h  =  r  sin  0  =  (r,-j  sin  0ij  +  rl>J+1  sin  0,j+i)  /2  (2.41) 

f  =  r  cos  0  =  cos  +  r,j+i  cos  0;  j+i)  /2  (2.42) 

and  at  the  Vj-j  node: 

(/  -  f)h  =  r  sin  0  =  (r,j  sin  0tJ  +  r.+ij  sin  0,+ij)  /2  (2.43) 

f  =  r  cos  0  =  ( ritj  cos  0,-j  +  r,+1j  cos  0t+i,i)  /2  (2.44) 


Nothing  has  been  said  up  to  this  point  about  choosing  the  frequency  with  which  the  drag 
tensor  field  is  updated  using  the  DRM.  However,  the  overall  efficiency  of  the  2|-D  model  clearly 
depends  on  the  choice  of  update  frequency.  In  fact,  comparison  tests  show  that  the  relative 
computational  time  between  the  DRM  and  DAM  is  roughly  equal  to  40/m  where  m  is  the 
number  of  time  steps  between  updates.  In  other  words,  a  run  of  the  2|-D  model  with  m  =  4 
spends  10  seconds  of  CPU  time  running  the  DRM  for  every  1  second  of  CPU  time  running  the 
DAM. 

A  primary  consideration  when  choosing  the  update  frequency  is  that  the  bottom  stress 
be  able  to  respond  appropriately  to  changes  in  forcing.  There  are  only  a  few  ways  in  which 
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Basin 

Period,  hr:min 

Lake  Ontario 

4:48 

Lake  Michigan 

9:03 

Lake  Superior 

8:00 

Lake  Erie 

14:18 

Narragansett  Bay 

0:44 

San  Francisco  Bay 

1:55 

Table  2.4:  Fundamental  seiching  periods  which  might  be  interpreted  as  t  pressure  for  2|-D. 
model  runs  in  the  listed  basin  [after  Wilson  (1972)]. 

information  about  changes  in  forcing  can  be  transmitted  to  the  sea  bottom.  The  mechanisms 
accounted  for  by  the  2  |-D  model  are  propagation  of  pressure  signals  in  the  DAM  and  turbulent 
diffusion  of  surface  stress  in  the  DRM.  Ideally,  the  update  frequency  should  be  chosen  to  resolve 
adequately  the  effect  of  these  processes  on  drag  tensor  variation. 

Time  scales  for  the  propagating  pressure  signals  in  the  DAM  are  intimately  dependent  on 
the  bottom  and  the  geometry  of  the  model  domain.  The  fastest  sea  surface  signals  travel  with 
the  shallow  water  gravity  wave  speed,  y/gK,  so  their  time  scale  can  be  thought  of  as  being 
determined  by  the  horizontal  (£)  and  vertical  (%)  scales  of  the  basin.  Explicitly, 

£ 

t pressure  oc  (2.45) 

Examples  of  £  might  be  the  shelf  width  for  coastal  ocean  flows  or  the  basin  diameter  for 
circulation  within  enclosed  or  semi-enclosed  bodies  of  water.  The  time  scale  defined  by  the 
latter  is  related  to  the  familiar  seiching  period  described  in  many  textbooks  (e.g.  Lamb,  1953). 
Wilson  (1972)  presents  a  variety  of  observed  fundamental  seiching  periods  in  well-known  basins. 
Table  2.4  summarizes  a  few  of  these.  Clearly,  basin  size  and  water  depth  greatly  affect  the 
seiching  period. 

Earlier,  it  was  stated  that  the  use  of  depth-independent  time  derivatives  in  the  DRM  con¬ 
strains  the  maximum  frequency  for  resolvable  processes.  Since  the  DRM  does  not  accelerate 
different  levels  of  the  water  column  at  different  rates,  it  can  only  resolve  processes  which  vary 
on  a  time  scale  longer  than  the  time  it  takes  shear  stress  related  signals  to  diffuse  throughout 
the  water  column.  In  other  words,  the  water  column,  as  modelled  by  the  DRM,  should  be 
considered  as  a  time-averager  for  processes  with  periods  shorter  than  the  turbulent  diffusion 
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u*s,cm/s 

Diffusion  time,  hr:min 

1 

■ 

2 

4 

| 

Table  2.5:  Diffusion  times  for  various  shear  velocities  in  20m  of  water 


time-scale.  One  estimate  of  this  time  scale  can  be  made  using  simple  scaling  arguments  (e.g. 
Csanady,  1984).  If  the  diffusion  rate  is  assumed  to  be  proportional  to  the  boundary  shear 
velocity,  v%  =  |r|/p  ,  the  time  for  a  signal  to  traverse  the  water  column  is  0 (ft/u*). 

The  ratio  of  depth  to  shear  velocity  is  a  convenient  but  crude  first  estimate.  Certainly  a 
model  which  updates  the  drag  tensor  every  5 h/u*  is  computationally  more  efficient  than  one 
which  updates  it  every  /i/u*.  Therefore,  a  more  specific  time  scale  is  sought.  Madsen’s  (1977) 
Ekman  layer  model  can  be  used.  The  model  yields  a  solution  for  the  spin-up  of  a  surface  Ekman 
layer  due  to  an  impulsively  started  constant  wind  stress  over  the  open  ocean.  The  model  has 
an  eddy  viscosity  that  increases  linearly  away  from  the  free  surface  with  the  same  origin  and 
slope  as  the  upper  function  of  (2.24).  The  resulting  time-  and  depth-dependent  velocity  profile 
is  given  by 

.  •  t*  (  \  cos/t  +  ism  ft  j7  ,n 

u  +  tv  =  —  /  exp  I - =  ) - = - at  (2.46) 

k  Jo  \  KU+at )  t 

where  z  is  the  distance  from  the  free  surface,  and  the  wind  stress  is  applied  in  the  positive 
^-direction. 

A  diffusion  time  scale  can  be  obtained  from  (2.46)  by  specifying  a  threshhold  level  for  the 
velocity  magnitude  at  z  =  h.  Reaching  the  threshhold  should  indicate  that  fluid  at  z  =  h  is 
“aware”  that  a  wind  stress  has  been  applied  at  z  =  0.  The  threshold  is  arbitrarily  chosen  here 
as  |u+iu|  =  «*,.  In  other  words,  it  is  assumed  that  the  water  column  responds  to  changes  in  the 
wind  stress  in  the  time  it  takes  the  velocity  at  z  =  h  to  spin  up  from  rest  to  a  magnitude  equal  to 
the  surface  shear  velocity.  Although  (2.46)  appears  complicated,  the  solution  for  |u  +  iv\  =  u»3 
is  well  approximated  by 


^  ~  ^diffusion  —  ^ 


u. 


(2.47) 


So  the  Ekman  layer  model  suggests  a  time  scale  four  times  longer  than  h/u»4.  Example  diffusion 
times  calculated  from  (2.47)  for  h  =  20m  are  given  in  Table  2.5. 
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Both  pressure  and  diffusion  time  scales  must  be  considered  when  choosing  an  update  fre¬ 
quency.  Comparison  of  tables  2.4  and  2.5  indicates  that  either  one  might  set  the  update  time 
scale.  Once  the  time  scale  is  set  the  number  of  time  steps  between  updates  is  given  by 


m  — 


_  ^update 
At 


(2.48) 


It  is  important  to  note  that  the  appropriateness  of  a  particular  choice  of  Update  *s  subject  to 
the  reasonability  of  the  2|-D  model  results  it  produces. 
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Chapter  3 


2^-D  Model  Results 


This  chapter  contains  2|-D  model  results  for  four  different  wind-driven  coastal  flow  scenarios. 
Each  has  a  domain  geometry  and  wind  stress  field  chosen  to  emphasize  a  particular  aspect  of 
the  relationship  between  depth-averaged  flow  and  bottom  stress.  The  domains  are  periodic  in 
the  alongshore  direction  and  bordered  by  a  vertical  drop-off  to  the  deep  ocean  at  a  straight 
offshore  boundary.  This  is  done  in  order  to  minimize  numerical  boundary  effects  which  might 
otherwise  contaminate  the  results.  The  wind  stress  field  is  also  chosen  to  be  spatially  uniform 
in  all  scenarios  for  the  same  reason. 

The  four  scenarios  are  (figure  3-1): 

1.  a)  Spin-up  from  rest  of  a  constant  depth  shelf  adjacent  to  a  straight  coast  for  a  steady 
wind  stress. 

2.  b)  Spin-up  from  rest  of  a  linearly-varying  depth  shelf  adjacent  to  a  straight  coast  for  a 
steady  wind  stress. 

3.  c)  Spin-up  from  rest  of  a  domain  with  sinusoidal  isobaths  adjacent  to  a  sinusoidal  coast 
for  a  steady  wind  stress. 

4.  d)  Response  of  a  linearly-varying  depth  shelf  adjacent  to  a  straight  coast  for  temporally 
periodic  wind  stress. 
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3.1  Constant  Depth  Shelf  Adjacent  to  a  Straight  Coastline 


Flow  over  a  flat  shelf  adjacent  to  a  long  straight  coastline  is  the  simplest  scenario  from  which 
relevant  conclusions  can  be  drawn.  These  results  are  used  both  to  introduce  the  physical 
phenomena  encountered  in  wind-driven  circulation  near  an  open  coast  and  to  examine  the 
basic  behavior  of  the  bottom  stress  -  depth-averaged  flow  relationship. 

3.1.1  Model  domain  and  wind  stress  field 

Specifically,  this  case  employs  a  rectangular  model  domain  (figure  3.1)  forced  by  a  wind  stress 
applied  uniformly  over  the  region  at  time  t  =  0  and  held  constant  thereafter.  The  bathymetry  is 
shallow  and  constant  ( h  =  20m)  between  y  =  0  and  y  =  —L.  A  no-flux  condition  perpendicular 
to  the  coastal  boundary  is  imposed  along  y  —  0  and  a  clamped  surface  height  condition  is 
imposed  along  y  =  —L  (as  if  the  domain  suddenly  became  infinitely  deep  there).  The  shelf’s 
cross-shore  dimension  is  195  km  while  its  alongshore  extent  is  70  km.  The  bottom  is  assumed  to 
be  covered  by  roughness  elements  with  corresponding  roughness  lengths,  zq  =  1  cm.  This  value 
is  chosen  so  as  to  represent  a  typical  roughness  length  due  to  wave-current  interaction.  Lastly, 
the  domain  is  overlain  by  an  Arakawa  “C”  grid  with  a  grid  length,  Ax  =  10  km.  Choosing  the 
grid  size  sets  the  CFL  time  step,  as  given  by  the  righthand  side  of  (2.21),  to  be  Atmax  =  505 
s  for  all  of  the  flat  bottom  runs.  The  actual  model  time  step,  however,  is  chosen  as  90%  of  the 
CFL  time  step,  or  455  s,  to  help  ensure  stability. 

The  short  alongshore  extent  of  the  domain  is  feasible  because  of  the  assumed  alongshore 
uniformity  in  the  response.  Since  periodic  boundary  conditions  are  implemented  at  the  cross- 
shore  boundaries,  end  effects  from  improperly  transmitted  waves  are  eliminated. 

While  a  wide  flat  shelf  is  clearly  a  drastic  simplification  of  realistic  shelf  topographies, 
this  first  scenario  illustrates  much  of  the  underlying  physics  of  wind-driven  flow  near  an  open 
coast  and  emphasizes  some  differences  between  the  2|-D  model  and  traditional  bottom  stress 
formulations.  More  realistic  bathymetries  and  wind  stress  fields  are  examined  later. 

Three  wind  stress  magnitudes  are  considered  in  this  first  scenario.  The  strongest  corresponds 
to  a  wind  speed  of  33  m/s  (63  knots)  at  10m  above  the  sea  surface.  This  yields  a  value 
of  ums/fh  =  30  (see  Table  2.1).  The  intermediate  strength  wind  stress  corresponds  to  a  wind 
speed  of  14  m/s  (27  knots)  yielding  a  value  of  uta/fh  =  10.  The  weakest  wind  stress  corresponds 
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to  a  wind  speed  of  8  m/s  (15  knots)  yielding  u*sf fh  =  5.  These  cases  will  be  referred  to  as 
the  shallow,  intermediate  and  deep  cases  respectively.  This  is  in  reference  to  the  relationship 
between  the  surface  boundary  layer  thickness  and  the  water  depth,  not  the  water  depth  itself 
which  is  constant. 

The  diffusion  time  scales,  as  determined  by  (2.47),  are  1250,  3750  and  7500s.  In  all  cases, 
the  diffusion  time  scale  is  shorter  than  the  pressure  time  scale,  L/y/gK  —  13900s.  Therefore, 
the  vertical  diffusion  times  set  the  update  frequencies.  The  update  frequencies  used  are  every 
5,  8  and  16  time  steps  for  the  shallow,  intermediate  and  deep  cases,  respectively. 

In  order  to  introduce  the  various  phenomena  efficiently,  a  case  with  both  alongshore  and 
cross-shore  wind  components  is  considered.  Until  specified  later,  all  of  the  results  presented 
are  for  a  wind  stress  directed  onshore  at  an  angle  of  45°  (measured  counter-clockwise  from  the 
x-axis). 


3.1.2  Surface  height  and  depth-averaged  velocity 

Figure  3-2  shows  time  series  of  sea  surface  height  for  the  three  wind  strengths.  Significant 
variation  among  the  subplots  illustrates  clearly  that  the  magnitude  of  the  wind  stress  has  a 
large  effect  on  sea  surface  response.  Notably,  there  are  differences  in  both  the  transient  and 
steady-state  behavior  of  the  flow. 

Roughly  speaking,  there  are  two  kinds  of  time-dependent  behavior  depicted  in  figure  3-2. 
The  first  has  a  long  smooth  variation  relative  to  the  second  and  represents  the  overall  approach 
to  steady-state.  The  second  has  a  shorter  more  oscillatory  variation  and  represents  the  wavelike 
adjustment  process  accompanying  the  spin-up.  Both  types  of  motions  are  present  in  the  three 
cases,  but  to  different  extents. 

The  spin-up  consists  of  a  short  initial  time  period  over  which  the  surface  remains  flat.  This 
lasts  at  each  point  until  the  arrival  of  the  signal  propagating  outward  from  the  coast.  The  signal 
travels  with  the  long  wave  phase  speed,  c  =  y/gK  =  14  m/s.  Therefore,  the  time  for  it  to  reach 
a  particular  point  is  a  function  of  that  point’s  distance  from  the  coast.  For  the  onshore  point 
(located  75  km,  or  7.5Ax,  from  the  coast),  the  time  should  be  75000  m  /  14  m/s  =  5360s  «  1.5 
hrs.  For  the  offshore  point  (located  135  km,  or  135Ax,  from  the  coast),  it  should  be  135000 
m  /  14  m/s  =  9640s  «  2.7  hrs.  It  is  difficult  to  check  this  in  figure  3-2.  However,  the  actual 
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Figure  3-2:  Time  series  of  surface  height  at  two  different  points  across  the  shelf.  The  point 
labelled  “inshore”  is  75  km  (7.5Ax)  from  the  shore  and  the  point  labelled  “offshore”  is  135  km 
(13.5Ax)  from  the  shore. 
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numerical  results  show  these  times  to  be  correct.  It  is  possible  to  see  that  the  onshore  point  in 
each  case  is  the  first  to  show  a  non-zero  height. 

Once  the  surface  height  becomes  non-zero,  all  three  cases  experience  an  initial  set-up  as 
the  cross-shore  component  of  wind  stress  moves  water  onshore.  Overall  rising  of  the  surface 
continues  only  in  the  shallow  case.  In  the  intermediate  and  deep  cases,  the  set-up  eventually 
stops  and  water  begins  to  move  offshore  producing  a  steady-state  set-down.  The  reason  for  this 
is  that,  at  steady-state,  the  Coriolis  term  predominates  in  the  cross- shore  momentum  balance 
for  the  intermediate  and  deep  cases  while  the  wind  stress  term  predominates  in  the  shallow 
case.  The  surface  slope  is  merely  a  response  that  balances  the  dominant  term  in  each.  The 
time  it  takes  for  the  overall  approach  to  steady-state  varies  from  about  5  hrs  in  the  shallow 
case  to  10  hrs  in  the  intermediate  case  to  20  hrs  in  the  deep  case,  which  is  an  indication  that 
frictional  damping  is  larger  for  larger  u+,/fh. 

Superimposed  upon  the  long  timescale  spin-up  are  the  higher  frequency  coastal  seiches. 
The  first  “crest”  of  these  waves  is  the  initial  set-up  propagating  out  from  the  coast.  The  signal 
reflects  off  the  clamped  offshore  boundary  and  travels  onshore  until  reflecting  off  the  coast. 
Reverberation  continues  until  the  seiches  are  damped  out.  The  time  series  in  figure  3-2  are 
records  of  these  damped  waves  as  they  pass  particular  points  on  the  shelf. 

The  plots  clearly  show  that  the  decay  timescale  for  the  waves  is  a  function  of  wind  strength. 
The  shallow  case  has  the  shortest  damping  timescale  and  the  deep  case  the  longest.  Therefore, 
the  damping  time  scales  axe  ordered  the  same  way  as  the  time  scales  for  the  overall  approach 
to  steady-state.  Both  relationships  are  an  indication  of  enhanced  frictional  effects  for  large 
um3/fh. 

It  is  easier  to  visualize  the  sea  surface  response  by  examining  figures  3-3-3-5  and  3-6-3-8 
instead  of  3-2.  Here,  surface  height  is  presented  as  a  function  of  cross-shore  distance  and  time. 
The  first  group  of  figures  displays  the  information  as  contour  lines  of  equal  surface  elevation 
and  the  second  as  a  three  dimensional  surface  formed  by  the  juxtaposition  of  “snapshots”  of 
cross-shore  sea  surface  profiles. 

Figures  3-3-3-5  show  very  clearly  the  initial  set-up  and  propagation  of  the  signal  away  from 
the  coast  at  early  times.  The  speed  of  propagation  can  roughly  be  deduced  by  taking  the  slope 
of  one  of  the  leftmost  contour  lines.  For  instance,  the  10  cm  contour  line  in  figure  3-3  has  a 
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Figure  3-3:  Contours  of  sea  surface  height  ii 
is  At  on  the  horizontal  axis  and  Ax  on  the 
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Figure  3-6:  3-D  plot  of  sea  surface  height  in  y  —  t  space  for  the  shallow  case.  Tic  mark  spacing 
is  At  on  the  time  axis  and  Ax  on  the  y  axis. 


Figure  3-7:  3-D  plot  of  sea  surface  height  in  y  —  t  space  for  the  intermediate  case.  Tic  mark 
spacing  is  At  on  the  time  axis  and  Ax  on  the  y  axis. 
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Figure  3-8:  3-D  plot  of  sea  surface  height  in  y  —  t  space  for  the  deep  case.  Tic  mark  spacing  is 
At  on  the  time  axis  and  Ax  on  the  y  axis. 

slope  of  about  190  km  /  30At  yielding  a  propagation  speed  of  13.9  m/s.  This  is  very  close  to 
the  theoretical  speed,  14.0  m/s.  Also  apparent  from  the  figures  is  the  reflection  of  the  initial 
sea  surface  signal  from  the  offshore  boundary  and  its  return  across  the  shelf  with  a  phase  speed 
equal  to  —y/gh.  The  number  and  relative  strength  of  the  reflections  vary  among  the  cases 
according  to  the  strength  of  frictional  damping. 

In  figures  3-4  and  3-5,  the  eventual  set-down  is  indicated  by  the  presence  of  negative  contours 
for  large  times.  Figure  3-3,  on  the  other  hand,  shows  only  positive  contours  indicating  that  the 
sea  surface  never  falls  below  its  initial  level.  The  value  of  the  steady-state  surface  slope  can 
be  roughly  deduced  from  the  figures  by  measuring  the  distance  between  contour  lines  at  The 
rightmost  edge  of  each  plot.  This  gives  a  steady-state  surface  slope  of  5.5  X  10-6  for  the  shallow 
run  and  set-down  surface  slopes  of  —1.8  X  10-6  and  —1.3  X  10-6  for  the  intermediate  and  deep 
runs.  The  latter  two  are  not  yet  at  steady-state,  but  longer  model  runs  show  that  their  surface 
slopes  are. 

An  equivalent  and  more  dramatic  way  of  viewing  the  evolution  of  the  free  surface  is  shown 
in  figures  3-6-3-8.  It  clearly  shows  both  the  steady-state  behavior  and  the  reverberations.  The 
latter  shows  up  as  the  characteristic  zig-zag  pattern  (also  seen  in  figures  3-3-3-5).  Even  though 
the  vertical  scale  in  figure  3-6  is  different  from  that  of  figures  3-7  and  3-8,  it  is  easy  to  tell  that 
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the  shallow  case  has  the  largest  response.  In  addition,  it  has  the  fastest  overall  spin-up  and  its 
reverbarations  die  out  the  fastest.  Both  of  these  are  evidence  of  the  larger  frictional  influence 
in  the  shallow  case. 

The  wavelike  structure  in  the  signal  is  more  obvious  in  figures  3-6-3-8  than  in  the  contour 
plots.  For  example,  in  the  shallow  case,  there  is  a  smaller  wave  propagating  out  of  phase  with 
the  one  originating  at  the  coast.  This  second  wave  originates  at  the  shelf  break  and  travels 
onshore  at  first.  It  is  not  visible  in  the  contour  plots,  but  is  seen  clearly  in  figure  3-6  as  a  second 
set  of  zig-zags.  The  wave  is  created  by  a  small  flux  divergence  due  to  application  of  the  clamped 
boundary  condition  at  the  shelf  break.  There  are  other  wavelike  structures  propagating  in  all 
of  the  runs.  They  are  responsible  for  some  of  the  additional  variability  between  the  onshore 
and  offshore  points  in  each  subplot  of  figure  3-2. 

Figure  3-9  shows  time  series  of  cross-shore  depth- averaged  velocity  for  the  three  different 
wind  stresses.  Like  figure  3-2,  there  is  a  significant  difference  in  transient  behavior  among  the 
subplots.  Unlike  figure  3-2,  however,  there  is  no  difference  in  steady-state.  In  all  cases,  the 
cross-shore  velocity  is  initially  onshore  as  the  water  moves  downwind.  It  remains  almost  entirely 
positive  in  the  shallow  case  as  the  sea  surface  sets  up.  The  intermediate  and  deep  cases  show 
a  large  offshore  flow  shortly  after  the  initial  set-up.  Oscillation  between  onshore  and  offshore 
flow  continues,  particularly  in  the  deep  case.  The  oscillations  are  damped,  however,  so  there  is 
a  net  offshore  flow  creating  a  set-down  in  each  case. 

All  of  the  runs  must  have  a  steady-state  cross-shore  velocity  equal  to  0  in  order  to  satisfy 
mass  conservation.  This  is  shown  in  (2.3)  which,  with  no  alongshore  variation,  reduces  to 
dV/dy  =  0  at  steady-state.  Since  V  is  always  0  at  the  coast,  it  must  be  0  everywhere. 

The  alongshore  depth-averaged  velocity  signal  is  shown  in  figure  3-10.  As  expected,  there 
is  significant  variability  among  the  cases.  The  relative  spin-up  times  and  oscillation  decay 
times  show  the  same  qualitative  behavior  as  the  surface  height  and  cross-shore  velocity.  At 
steady-state,  the  shallow  case  predictably  has  the  largest  alongshore  velocity,  followed  by  the 
intermediate  and  deep  cases. 

Taken  together,  the  sea  surface  and  depth-averaged  velocity  behavior  indicates  some  fun¬ 
damental  properties  which  have  implications  for  this  thesis.  They  are: 

•  Many  of  the  differences  among  the  cases  are  due  to  differences  in  the  relative  importance 
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Figure  3-9:  Time  series  of  cross-shore  depth-averaged  veloctiy  at  two  different  points  across 
the  shelf.  The  point  labelled  “inshore”  is  75  km  (7.5Ax)  from  the  shore  and  the  point  labelled 
“offshore”  is  135  km  (13.5 Ai)  from  the  shore. 
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Figure  3-10:  Time  series  of  alongshore  depth-averaged  veloctiy  at  two  different  points  across 
the  shelf.  The  point  labelled  “inshore”  is  75  km  (7.5Ax)  from  the  shore  and  the  point  labelled 
“offshore”  is  135  km  (13.5 Ax)  from  the  shore. 
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of  bottom  stress.  This  is  supported  by  the  fact  that  the  differences  in  transient  behavior 
primarily  consist  of  variations  in  spin-up  time  and  damping  of  the  oscillations,  both  of 
which  are  frictional  in  nature. 

•  The  magnitude  of  u*s//h  is  a  significant  factor  in  determining  the  importance  of  bottom 
stress  in  wind-driven  dynamics  near  an  open  coast  because  the  only  difference  among  the 
cases  is  their  assumed  values  of  u*s/fh.  The  results  of  Jenter  and  Madsen  (1989)  also 
support  this  idea. 

3.1.3  Momentum  balances 

The  assertions  are  confirmed  by  further  examining  the  dynamics  of  the  shallow,  intermediate 
and  deep  cases.  This  subsection  contains  detailed  descriptions  of  the  momentum  balances  for 
each.  Only  results  for  the  offshore  point  are  presented  because  there  is  no  qualitative  difference 
among  locations  across  the  shelf. 

In  all  three  cases  (figures  3-lla,b;  3-12a,b  and  3-13a,b),  the  very  early  momentum  balance  is 
for  a  forced  inertial  oscillation.  The  first  time  step  requires  the  local  time  derivative  to  balance 
the  wind  stress  (since  the  Coriolis  term,  pressure  gradient  and  bottom  stress  are  all  0  at  t  =  0). 
As  a  downwind  flow  develops,  the  Coriolis  term  enters  the  equation.  Pure  inertial  oscillations 
last  only  until  the  bottom  stress  first  gets  updated.  This  happens  at  t  «  40  min,  f  »  1  hr  and 
t  «  2  hrs  in  the  shallow,  intermediate  and  deep  cases,  respectively. 

As  the  bottom  stress  develops,  the  degree  to  which  the  inertial  oscillations  are  affected 
varies.  Frequencies  and  decay  scales  become  a  function  of  shelf  geometry  and  strength  of 
bottom  friction  (e.g.  Madsen,  in  preparation).  In  the  shallow  case,  the  oscillations  are  strongly 
damped  (figures  3-lla,b)  while  in  the  deep  case  they  persist  for  a  relatively  long  time  (figures 
3.13a, b),  but  with  a  frequency  shift. 

In  all  three  cases,  there  is  a  jagged  quality  to  the  alongshore  bottom  stress  and  acceleration 
signals.  Jumps  occur  each  time  the  bottom  stress  is  updated,  but  are  more  pronounced  for 
smaller  u*„/fh.  Particularly  in  the  intermediate  case,  the  jumps  seem  to  straddle  the  basic 
signal,  alternating  each  time  the  bottom  stress  is  updated.  This  is  not  an  indication  of  a  poorly 
resolved  physical  process,  but  rather  the  expression  of  small  numerical  waves  produced  and 
reinforced  each  time  the  bottom  stress  is  updated.  These  waves  also  appear  in  the  cross-shore 


54 


cm/s‘»2  cm/s 


time,  hrs 


time,  hrs 


Figure  3-11:  Time  series  of  shallow  case  momentum  balance  components  at  the  offshore  point. 
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Figure  3-12:  Time  series  of  intermediate  case  momentum  balance  components  at  the  offshore 
point. 
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Figure  3-13:  Time  series  of  deep  case  momentum  balance  components  at  the  offshore  point. 
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momentum  balance  in  the  pressure  gradient  and  local  time  derivative  terms.  Clearly,  the  choice 
of  update  period  has  a  strong  influence  on  their  nature.  This  emphasizes  the  fact  that  the  update 
time  scale  must  be  separate  from,  and  much  shorter  than,  the  physical  time  scales  of  interest 
in  order  for  the  2^-D  model  to  be  appropriate.  Care  must  be  taken  in  each  application  of  the 
model  to  determine  the  time  scales  of  interest  and  to  choose  the  update  frequency  accordingly, 
withing  the  constraints  created  by  approximating  the  local  time  derivatives  in  the  DRM. 

The  alongshore  component  of  bottom  stress  in  all  cases  grows  until  it  completely  balances  the 
alongshore  wind  stress  (figures  3-lla,  3-12a  and  3-13a).  This  means  the  steady-state  alongshore 
flow  is  in  Couette  balance.  It  is  easy  to  understand  this  by  examining  (2.1).  The  only  terms 
left  at  steady-state  are  the  bottom  and  surface  stresses.  Because  of  the  uniformity  assumed 
here,  the  alongshore  pressure  gradient  is  always  zero.  Also,  the  Coriolis  term  must  go  to  zero 
as  V  goes  to  zero.  Thus,  the  alongshore  bottom  stress  plays  a  major  role  in  determining  the 
alongshore  steady-state  flow. 

The  relative  importance  of  cross-shore  bottom  stress  in  determining  the  cross-shore  steady- 
state  flow  is  variable.  Unlike  the  alongshore  bottom  stress,  it  is  not  constrained  by  the  lack  of 
other  forces  to  balance  the  wind  stress.  Therefore,  its  role  is  dependent  on  the  strength  of  the 
other  terms.  From  figures  3-llb,  3- 12b  and  3-13b,  it  is  clear  that  the  Coriolis  term,  surface 
stress,  pressure  gradient  and  bottom  stress  all  play  a  role  in  the  balance. 

In  the  shallow  case,  the  Coriolis  term  is  slightly  less  than  the  surface  stress  and  directed 
offshore.  The  pressure  gradient  is  directed  offshore  also,  corresponding  to  the  steady-state  set¬ 
up  (e.g.  see  figure  3.1).  The  bottom  stress,  though  not  a  major  term  in  the  momentum  balance 
is  finite  and  in  the  same  direction  as  the  wind  stress.  Thus  all  of  the  terms  contribute  to  the 
balance. 

It  is  worthwile  to  repeat  at  this  time  that  a  non-zero  steady-state  cross-shore  bottom  stress 
cannot  be  predicted  by  a  DAM  using  a  traditional  bottom  stress  formulation.  Therefore,  the 
momentum  deficit  created  by  ignoring  the  cross-shore  bottom  stress  would  be  made  up  for  in  a 
traditional  DAM  by  changing  the  other  terms.  Since  the  alongshore  flow  is  predetermined  by 
the  Couette  balance,  the  compensation  wouldt  be  made  by  altering  the  pressure  gradient.  This 
implies  that  the  neglect  of  cross-shore  bottom  stress  directly  affects  the  prediction  of  cross-shore 
surface  slope  (and  its  integral,  surface  height).  A  quantification  of  these  differences  is  presented 
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later. 

The  Coriolis  term  and  pressure  gradient  play  a  progressively  more  important  role  in  the 
cross-shore  momentum  balance  as  the  wind  gets  weaker.  Conversely,  bottom  stress  plays  a 
decreasingly  important  role.  This  indicates  merely  that  the  cross-shore  balance  is  becoming 
more  geostrophic  as  the  wind  stress  gets  weaker. 

As  the  bottom  stress  decreases  in  relative  importance  with  lessening  wind  stress,  its  neglect 
in  the  cross-shore  direction  requires  less  adjustment  of  the  cross-shore  pressure  gradient.  Con¬ 
sequently,  there  is  little  difference  in  sea  surface  height  signals  between  the  2^-D  model  and  a 
traditional  DAM  for  very  weak  winds. 

3.1.4  Bottom  stress  and  drag  tensor 

Figures  3-14-3-16  show  time  series  of  bottom  stress  for  the  different  values  of  u*s/fh.  The 
continuous  lines  are  from  the  2  |-D  model1  and  the  symbols  are  estimates  produced  by  the 
DRM.  The  horizontal  spacing  of  the  symbols  and  the  jaggedness  of  the  continuous  curves  in 
each  plot  indicate  the  update  period  for  the  drag  tensor.  The  reason  for  the  jaggedness  is 
that  between  updates  the  drag  tensor  is  held  constant.  In  reality,  both  the  drag  tensor  and 
depth-averaged  velocity  should  change  continuously.  However,  since  the  tensor  is  constant 
between  updates,  only  the  depth-averaged  velocity  alters  the  bottom  stress.  This  is  insufficient 
to  account  completely  for  the  true  bottom  stress  variation.  Therefore,  at  the  next  update  time, 
the  drag  tensor  makes  an  abrupt  change  and  the  DAM  bottom  stress  jumps  accordingly,  as  was 
seen  earlier. 

Despite  the  jaggedness,  the  symbols  and  continuous  lines  agree  very  well.  This  is  important 
because  they  are  not  entirely  constrained  to  do  so.  Agreement  is  particularly  good  near  steady- 
state.  The  better  agreement  at  steady-state  is  likely  an  indication  that  treatment  of  time 
dependence  in  the  DRM  is  not  completely  adequate,  which  is  not  surprising  since  the  DRM 
was  simplified  by  assuming  depth-independent  time  derivatives. 

Agreement  is  better  in  the  shallow  case  than  in  the  other  two.  This  likely  indicates  that  the 
barotropic  time  derivatives  are  a  better  approximation  in  shallow  water,  and  as  discussed  above, 

1Note  that  the  continuous  lines  are  merely  the  bottom  stress  curves  from  figures  3-11-3-13  multiplied  by  the 
water  depth. 
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Figure  3-14:  Time  series  of  shallow  case  kinematic  bottom  stress  components  Tbx/p,Tly/p  from 
the  DAM  (continuous  curves)  and  from  the  DRM  (symbols). 
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BOTTOM  STRESS  -  DAM  and  DRM,  INTERMEDIATE 


Figure  3-15:  Time  series  of  intermediate  case  kinematic  bottom  stress  components  Tbx/p,Tby/p 
from  the  DAM  (continuous  curves)  and  from  the  DRM  (symbols). 
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Figure  3-16:  Time  series  of  deep  case  kinematic  bottom  stress  components  Ux/p^y/p  from  the 
DAM  (continuous  curves)  and  from  the  DRM  (symbols). 


that  the  length  of  the  update  period  affects  the  agreement.  The  two  are  clearly  related  since 
the  approximate  treatment  of  time  dependence  in  the  DRM  prohibits  arbitrarily  shortening  the 
update  period. 

Figure  3-17  shows  time  series  of  drag  tensor  magnitude.  It  is  apparent  that  the  temporal 
behavior  of  r  is  strongly  dependent  on  the  strength  of  the  wind  stress.  In  each  case,  r  starts 
out  at  zero  and  grows  to  its  steady-state  value  over  roughly  the  spin-up  time  mentioned  earlier. 
Both  the  shallow  and  intermediate  cases  experience  an  overshoot  in  the  resistance  coefficient 
before  reaching  steady-state.  The  steady-state  r  values  are  clearly  a  function  of  u*s/ fh  with 
the  shallow  case  being  largest. 

The  basic  time-varying  behavior  is  not  surprising,  if  one  believes  that  the  bottom  stress 
magnitude  is  better  approximated  by  a  quadratic,  rather  than  a  linear,  drag  law.  In  a  quadratic 
drag  law,  the  quantity  with  the  same  dimensions  as  the  resistance  coefficient  is  the  product  of 
the  drag  coefficient  and  the  depth-averaged  speed.  Therefore,  if  r  varies  linearly  with  the 
magnitude  of  the  depth-averaged  velocity,  it  is  merely  mimicking  a  quadratic  drag  law. 

This  can  be  tested  by  deriving  an  estimated  drag  coefficient  from  the  2  |-D  model  results. 
Figure  3-18  shows  the  ratio  of  the  resistance  coefficient  to  the  magnitude  of  the  depth-averaged 
velocity.  Interestingly,  the  curves  show  as  much  temporal  variability  as  their  counterparts 
in  figure  3-17.  However,  most  of  the  steady-state  variability  with  ums/fh  has  been  removed. 
Jenter  and  Madsen  (1989)  have  addressed  the -steady-state  problem  and  similarly  found  that 
the  steady-state  drag  coefficient  is  relatively  independent  of  u„s/ fh.  They  also  found  that  the 
main  parameter  determining  the  value  of  Cd  is  Zo/h.  Since  the  roughness  length  to  water  depth 
ratio  is  the  same  in  all  three  cases  here,  the  steady-state  drag  coefficient  is  roughly  the  same. 

It  is  somewhat  surprising  that  the  initial  overshoots  seen  in  figure  3-17  are  not  absent  in  the 
estimated  drag  coefficient  curves  since  there  are  corresponding  overshoots  in  the  velocity  signals 
(figures  3-9  and  3-10).  In  fact,  the  relative  overshoot  seems  to  be  enhanced  for  the  intermediate 
and  deep  cases  and  unchanged  for  the  shallow  case.  Apparently  in  the  region  of  the  overshoot, 
the  bottom  stress  magnitude  is  increasing  faster  than  the  depth-averaged  speed  squared.  Also 
surprising  is  the  jaggedness  of  the  intermediate  and  deep  case  Cd  signals.  Apparently,  the  small 
numerically-generated  waves  superimposed  on  the  physical  signal  are  capable  of  modulating  the 
drag  coefficient.  As  emphasized  throughout  this  thesis,  however,  only  processes  which  occur  on 
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Figure  3-18:  Time  series  of  estimated  drag  coefficient  at  the  offshore  point. 
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time  scales  much  longer  than  the  update  frequency  are  resolved  appropriately.  The  very  high 
frequency  variation  is  non-physical. 

The  estimated  drag  coefficient  results  imply  that  the  steady-state  bottom  stress  magnitude 
can  be  determined  fairly  accurately  by  employing  a  quadratic  drag  law  with  a  drag  coefficient 
based  on  a  good  estimate  of  Zo/h.  Unfortunately,  a  correct  estimate  of  bottom  stress  magnitude 
is  not  a  correct  estimate  of  bottom  stress.  The  direction  of  the  stress  still  needs  to  be  determined. 

Figures  3-19-3-21  depict  the  angular  information  necessary  to  determine  the  rotational  part 
of  the  drag  tensor  for  each  of  the  three  wind  strengths  at  early  times.  The  curve  labeled  “theta” 
(0)  is  the  angular  difference  between  the  bottom  stress  and  depth-averaged  flow  predicted  by  the 
DRM.  The  curve  labeled  “alpha”(a)  is  the  direction  of  the  bottom  stress  and  the  one  labeled 
“  beta”(/3)  is  the  direction  of  the  depth-averaged  flow.  These  angles  are  defined  in  figure  2-1. 

In  each  figure,  the  (3  curve  starts  out  as  expected  for  an  inertial  oscillation.  The  direction 
of  flow  is  linearly  proportional  to  t  until  the  bottom  stress  is  turned  on.  Rotation  continues 
relatively  unperturbed,  because  bottom  stress  is  initially  very  small,  until  the  cross-shore  pres¬ 
sure  gradient  arrives.  At  this  time,  the  frictionally-  and  geometrically-modulated  oscillations 
commence  and  the  angular  variation  deviates  from  linear.  The  directional  oscillations  persist 
for  a  time  commensurate  with  the  level  of  frictional  influence  in  each  case,  but  eventually  die 
down  such  that  the  steady-state  direction  becomes  0°  or  alongshore. 

The  9  curve  is  more  difficult  to  interpret,  particularly  during  the  first  few  update  periods. 
In  the  shallow  case,  0  starts  at  approximately  —75°,  but  flips  up  to  10°  at  the  second  update 
time.  This  constitutes  a  change  of  roughly  80°  in  40  min.  This  is  not  likely  physical,  but  rather 
a  problem  induced  by  the  DRM.  In  fact,  it  is  likely  a  matter  of  the  initial  stresses  being  so 
small  that  their  components  are  poorly  constrained  by  the  convergence  criteria  of  the  DRM. 

After  the  initial  jump,  however,  the  9  curves  behave  fairly  regularly.  In  the  shallow  and 
intermediate  cases,  the  bottom  stress  is  at  first  briefly  rotated  clockwise  from  the  depth-averaged 
flow.  However,  it  soon  passes  through  the  direction  of  the  flow  and  remains  rotated  slightly 
counter-clockwise  until  steady-state  is  reached.  This  means  that  since  the  steady-state  flow  is 
in  the  positive  x  direction,  the  bottom  stress  is  directed  slightly  offshore  as  was  indicated  in 
the  momentum  balances.  The  deep  case  bottom  stress  is  always  directed  clockwise  from  the 
depth-averaged  flow. 
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The  a  curve  is  just  the  sum  of  the  0  and  /3  curves.  In  other  words,  the  direction  of  the 
bottom  stress  in  the  2|-D  model  is  rotated  from  the  direction  of  the  depth-averaged  flow  by  a 
value  equal  to  the  instantaneous  value  of  0.  Since  0  is  only  updated  periodically,  the  a  and  (3 
curves  are  always  the  same  distance  apart  between  updates.  However,  the  distance  varies  each 
time  0  is  updated. 

Notice  that  figures  3-20  and  3-21  do  not  show  steady-state.  Only  early  times  are  displayed  in 
order  to  clearly  illustrate  the  initial  angular  behavior.  The  angular  information  reaches  steady- 
state  in  times  comparable  to  the  spin-up  times  for  each  case.  The  steady-state  value  of  (3  is  0 
for  all  the  cases  since  V  =  0.  The  steady-state  values  of  0  are  —6°,  —4°,  —1°  for  the  shallow, 
intermediate  and  deep  cases,  respectively.  The  progressive  decrease  in  0  is  an  indication  that 
the  cross-shore  bottom  stress  decreases  in  relative  importance  from  the  shallow  to  deep  cases 
since  0  =  a  =  arctan(r(,y/r(,x).  Figures  3-11-3-13  support  this. 

Another  thing  that  the  angular  results  illustrate  is  that  angular  differences  between  the 
bottom  stress  and  depth-averaged  flow  are  relatively  small.  Therefore,  it  seems  appropriate  to 
ask  “Does  including  them  in  a  drag  tensor  really  matter?”.  This  question  is  addressed  in  the 
next  subsection. 

3.1.5  Differences  between  traditional  and  2|-D  model  results 

For  an  initial  comparison,  a  traditional  model  was  run  using  a  resistance  coefficient  equal  to  the 
shallow  case  steady-state  resistance  coefficient  determined  from  figure  3-17.  Figure  3-22  depicts 
the  surface  height  variation  at  the  offshore  point  for  this  traditional  run  (dashed  line)  and  for 
the  2^-D  model  (solid  line).  There  is  clearly  a  difference  in  the  two  signals.  The  traditional 
model  spins  up  more  slowly  than  the  2|-D  model  which  makes  sense  in  light  of  the  actual 
time-dependent  behavior  of  r  shown  in  figure  3-17  (reproduced  as  a  subplot  in  figure  3-22). 
The  traditional  model  value  of  r,  essentially  equivalent  to  Cd|i7|steady-state>  produces  friction 
that  is  too  large  until  \U\  >  |^|steady-state-  The  larger  friction  also  accounts  for  increased 
smoothness  in  the  traditional  model  results. 

A  third  line  (dotted)  is  also  shown  in  figure  3-22.  It  corresponds  to  a  traditional  model 
using  a  quadratic  drag  law  with  Cd  =  ^steady-state  firom  the  2!'D  model.  The  constant  Cd 
model  spins  up  more  slowly  than  the  2|-D  model  and  more  rapidly  than  the  constant  r  model, 
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Figure  3-22:  Time  series  of  sea  surface  height  at  the  offshore  point  for  different  drag  law 
assumptions:  variable  drag  tensor  (solid  line),  constant  resistance  coefficient  (dashed  line), 
constant  drag  coefficient  (dotted  line). 
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Figure  3-23:  Time  series  of  alongshore  depth-averaged  velocity  at  the  offshore  point  for  different 
drag  law  assumptions:  variable  drag  tensor  (solid  line),  constant  resistance  coefficient  (dashed 
line),  constant  drag  coefficient  (dotted  line). 


commensurate  with  the  strength  of  bottom  friction  at  early  times. 

One  of  the  most  important  result  of  this  section  is  illustrated  by  the  difference  in  steady- 
state  surface  heights.  Because  the  steady-state  cross-shore  pressure  gradient  is  independent  of 
y  and  the  surface  is  clamped  at  y  —  —  L,  the  difference  between  traditional  and  2^-D  model 
surface  heights  at  any  point  is  linearly  proportional  to  the  distance  from  the  shelf  edge.  This 
means  the  ratio  of  traditional  and  2|-D  model  surface  heights  at  any  point  is  constant  across 
the  shelf.  Thus,  the  ratio  of  heights  at  the  offshore  point  can  also  be  interpreted  as  the  ratio  of 
coastal  set-ups.  The  difference  from  figure  3-22  is  19%. 

The  ability  of  the  small  rotation  of  the  bottom  stress  (—6°  in  the  shallow  case)  to  cause 
large  changes  in  sea  surface  slope  may  be  somewhat  surprising.  However,  as  discussed  in  the 
momentum  balance  subsection,  the  critical  question  is  not  “How  big  a  role  does  the  bottom 
stress  play  in  the  cross-shore  momentum  balance?”  but  “How  much  of  the  difference  between  the 
wind  stress  and  Coriolis  term  must  the  bottom  stress  support  relative  to  the  pressure  gradient?” 
As  can  be  seen  from  figure  3-llb,  the  steady-state  bottom  stress  represents  a  significant  fraction 
(roughly  20%)  of  the  steady-state  pressure  gradient. 

The  alongshore  depth-averaged  velocity  (figure  3-23)  spins  up  similarly  to  the  surface  height. 
The  traditional  models  are  retarded  because  of  their  initially  large  bottom  friction.  However, 
unlike  the  surface  height,  the  steady-state  values  are  almost  identical.  The  traditional  models’ 
velocities  are  only  slightly  smaller. 

The  agreement  between  steady-states  makes  sense  in  view  of  the  smallness  of  6.  For  example, 
the  traditional  model  resistance  coefficient  was  chosen  to  equal  the  steady-state  value  of  r  from 
the  2|-D  model.  Therefore,  the  alongshore  velocity  is  equal  to  tix /pr  cos  0  for  the  2|-D  model 
and  Tbx/pr  for  the  traditional  model.  Since  cos  0  «  1,  the  two  are  almost  the  same. 

The  alongshore  velocities  explain  why  the  difference  in  steady-state  pressure  gradient  be¬ 
tween  the  traditional  and  2|-D  models  is  nearly  equal  to  the  2|-D  model  cross-shore  bottom 
stress.  The  wind  stress  is  the  same  for  the  models,  as  is  the  Coriolis  term.  Therefore,  the 
momentum  deficit  to  be  made  up  by  the  pressure  gradient  and  bottom  stress  is  the  same. 
Whatever  is  not  balanced  by  bottom  stress  must  be  balanced  by  the  pressure  gradient. 

Figure  3-24  shows  the  cross-shore  velocity  for  the  traditional  and  2  |-D  model  runs.  All 
go  to  zero  at  steady-state,  as  they  must.  The  transient  response  mimicks  those  of  the  surface 
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Wind  Direction 

r,  cm/s 

9 

hahi 

0° 

.348 

■H 

.02 

30° 

.325 

BS 

-.78 

45° 

.298 

-6° 

-.19 

60° 

.256 

-9° 

-.13 

O 

o 

CM 

.261 

11° 

-.07 

135° 

.301 

-.07 

150° 

.328 

-.06 

180° 

.348 

-.02 

Table  3.1:  Differences  between  traditional  and  2|-D  model  steady-state  results  for  the  shallow 
case. 

height  and  alongshore  velocity. 

It  is  useful  to  examine  the  difference  in  surface  heights  between  the  models  more  closely. 
Figures  3-22-3-24  display  results  for  only  one  wind  direction.  Table  3.1  gives  ratios  of  steady- 
state  surface  slopes  for  various  wind  angles.  Notice  that  the  difference  between  the  models  is 
strongly  dependent  upon  wind  angle  but  the  order  of  magnitude  seems  to  be  roughly  0.1.  The 
very  large  underprediction  for  winds  at  30°  is  a  result  of  the  sea  surface  slope  being  almost 
zero  there.  Consequently,  the  ratio  is  large  because  the  denominator  is  small,  not  because  the 
numerator  is  large. 

There  is  an  asymmetry  between  shallow  cases  with  angles  less  than  90°  and  those  with 
angles  greater  than  90°.  Obtuse  wind  angles  show  smaller  differences  than  their  compliments 
because  Coriolis  accelerations  acts  with  the  wind  for  angles  greater  than  90°  requiring  larger 
set-ups.  The  differences  are  not  due  to  differences  in  bottom  stresses  because  bottom  stresses 
have  essentially  the  same  magnitude  and  direction  for  complimentary  wind  angles.  Therefore, 
increased  denominators  in  the  ratio  are  responsible  for  the  asymmetry  rather  than  increased 
absolute  differences. 

Notice  that  the  bottom  stress  is  almost  always  directed  offshore  from  the  depth-averaged 
velocity  when  the  wind  is  between  0°  and  180°.  This  has  important  implications  for  near  bottom 
transport  prediction.  It  means  that,  whenever  the  wind  has  an  onshore  component,  the  bottom 
stress  has  an  offshore  component.  Traditional  DAMs  can  predict  only  alongshore  steady-state 
transport.  Clearly,  9  cannot  be  ignored  if  near  bottom  transport  is  of  interest. 
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Figure  3-24:  Time  series  of  cross-shore  depth-averaged  velocity  at  the  offshore  point  for  different 
drag  law  assumptions:  variable  drag  tensor  (solid  line),  constant  resistance  coefficient  (dashed 
line)}  constat  dr«^coePP7c.ie/vt  faotlcd  Kr\e>. 
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Results  for  angles  near  90°  are  conspicuously  absent  from  table  3.1  because  the  alongshore 
momentum  balance  in  the  2  |-D  models  is  degenerate  when  the  wind  is  directly  onshore,  and 
alongshore  uniformity  of  the  domain  and  wind  stress  means  that  there  can  be  no  steady-state 
bottom  stress  component  in  the  alongshore  direction.  Since  there  can  be  no  depth-averaged 
cross-shore  flow,  9  must  be  90°.  The  alongshore  velocity  can  be  anything  when  this  is  the  case. 
Without  constraint  on  the  alongshore  flow,  it  can  easily  become  unstable.  In  numerical  practice, 
wind-stress  angles  between  85°  and  95°  tend  to  have  unstable  alongshore  flows.  The  instability 
is,  however,  not  a  problem  in  more  realistic  cases.  If  the  domain  geometry  or  wind  stress  has  any 
variability  in  the  cross-wind  direction,  a  pressure  gradient  will  be  set-up  to  balance  a  cross-wind 
bottom  stress  component. 

3.2  Linearly- Varying  Depth,  Straight  Coastline 

Flow  over  a  linearly-varying  depth  shelf  adjacent  to  a  long  straight  coastline  is  described  here. 
This  represents  an  increase  in  complexity  from  the  previous  case  where  there  was  no  cross¬ 
shore  depth  vaxiation.  Consequently,  the  two  dimensionless  parameters  which  determined  the 
importance  of  bottom  friction  in  the  flat  bottom  case,  u„s/ fh  and  zo/h,  are  now  variable  in  the 
cross-shore  direction.  Variation  implies  that  for  a  single  wind  stress,  this  scenario  combines  a 
range  of  behaviors  from  the  previous  section. 

3.2.1  Model  domain  and  wind  stress  field 

The  model  domain  for  this  scenario  is  rectangular  with  the  same  horizontal  dimensions  as  those 
in  section  3.1.  However,  the  bathymetry  now  varies  linearly  across  70%  of  the  basin  while  the 
outer  30%  is  flat.  The  depth  at  half  a  grid  space  from  the  coast  is  5  m  and  increases  by  5  m 
per  grid  space  (slope  =  5  x  10~4  since  Ax  is  again  10  km)  until  reaching  70  m.  The  outer  shelf 
depth  is  constant  at  70  m.  The  strip  of  flat  bathymetry  is  added  offshore  in  order  to  have  a 
depth  variation  similar  to  one  used  later  in  section  3.3.  It  is  important  to  note  that  there  are 
no  model  limitations  which  say  that  the  depth  profile  cannot  be  truncated  at  the  end  of  its 
linear  variation. 

The  bottom  is  again  assumed  to  be  covered  by  roughness  elements  corresponding  to  zq  =  1 
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cm  which  means  that  zo/h  varies  from  2  x  10-3  near  the  coast  to  1.4  x  10-4  offshore.  It  was 
5  x  10-4  in  the  flat  bottom  runs.  The  wind  stress  is  the  same  as  in  the  flat  bottom  shallow 
case.  It  corresponds  to  a  33  m/s  wind  speed  yielding  a  range  of  u*3/fh  from  8.6  in  deep  water 
to  120  at  the  coast.  The  wind  direction  is  45°.  Because  the  depth  is  variable,  the  CFL  time 
step  is  not  strictly  applicable.  However,  a  similar  criterion  is  applied.  The  maximum  allowable 
time  step  is  assumed  to  be  Afmax  =  Ax/V^/imax  =  270  s.  The  actual  time-step  used  is  90% 
of  this  value,  or  243  s. 

The  diffusion  time  scales,  as  determined  by  (2.47),  range  between  315  s  and  4375  s.  The  up¬ 
date  frequency  chosen  is  once  every  5  time  steps.  This  means  the  diffusion  process  is  somewhat 
under-resolved  near  the  coast  (depth  shallower  than  20m)  and  over-resolved  offshore  (depth 
deeper  than  20  m). 

3.2.2  Surface  height  and  depth-averaged  velocity 

Figure  3-25  shows  time  series  of  the  basic  variables  at  two  cross-shore  points.  The  onshore 
(offshore)  point  corresponds  to  a  water  depth  such  that  uma/ fh  equals  the  shallow  (intermediate) 
values  from  section  3.1.  Specifically,  the  onshore  point  is  in  20m  of  water  and  the  offshore  point  is 
in  60m  of  water.  Notice  that  although  u,4///i  at  the  offshore  point  agrees  with  the  intermediate 
case  of  section  3.1,  the  value  of  zo/h  is  smaller. 

Figure  3-25a  shows  that  the  surface  height  signal,  both  onshore  and  offshore,  looks  like  a 
combination  of  the  flat  bottom  shallow  and  intermediate  responses.  The  signal  magnitude  is 
comparable  to  the  flat  bottom  shallow  case,  but  the  spin-up  time  is  more  like  the  flat  bottom 
intermediate  case.  Reasons  for  this  are  given  below.  The  initial  response  at  both  points  is  flat, 
followed  by  a  time-lagged  set-up.  The  onshore  point  sets-up  first  reaching  a  value  comparable 
to  the  shallow  case  in  section  3.1.  This  happens  slightly  later,  however,  than  the  comparable 
point  in  the  flat  bottom  case  because  the  wave  travelling  out  from  the  coast  is  travelling  through 
shallower  water  here.  Conversely,  the  offshore  point  experiences  a  slightly  earlier  arrival  since 
the  increased  phase  speed  offshore  of  the  20m  depth  contour  allows  the  signal  to  more  than 
regain  the  time  lost  travelling  in  shallow  water.  Note  that  the  signal  at  the  offshore  point  has 
a  much  larger  amplitude  than  that  of  the  intermediate  case  in  section  3.1.  This  indicates  that 
the  amplitude  of  the  transient  wind-driven  signal  is  determined  at  the  coast  which,  makes  sense 
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because  the  initial  bulge  is  produced  there. 

After  the  short  set-up  period,  the  sea  surface  height  at  both  points  drops  significantly  as 
the  water  moves  olfshore.  This  creates  an  overall  set-down.  Damped  oscillations  continue 
superimposed  on  the  set-down  for  a  time  comparable  to  the  flat  bottom  intermediate  case 
decay  time,  whichis  understandable  because  the  waves  generated  at  the  coast  are  not  damped 
as  quickly  when  they  propagate  through  deep  water. 

Figure  3-26  presents  another  picture  of  the  spin-up  process.  The  wave  progresses  outward 
from  the  coast  in  much  the  same  way  as  the  flat  bottom  case,  but  the  spatially-varying  phase 
speed  causes  curvature  in  the  contour  lines.  The  phase  speed  is  slower  nearshore  and  faster 
offshore  giving  the  signature  a  cusp-like  appearance.  It  is  possible  to  determine  an  average 
phase  speed  based  on  the  total  travel  time  for  the  wave  to  cross  the  shelf  using  travel  time  = 
/^(l/cjdy.  For  the  present  domain  geometry,  this  yields  a  phase  speed  of  15.25  m/s.  A 
phase  speed  estimated  using  the  leftmost  10  cm  surface  contour  yields  an  average  phase  speed 
somewhere  between  14  and  16  m/s  which  is  in  good  agreement  with  the  theoretical  value. 

The  oscillations  eventually  damp  out,  and  a  good  approximation  to  steady-state  is  reached 
in  roughly  20  hrs.  Unlike  the  flat  bottom  cases,  the  steady-state  surface  slope  is  variable  across 
the  shelf.  It  changes  monotonically  from  negative  offshore  to  positive  onshore  with  the  steepest 
slope  found  at  the  coast.  Despite  the  positive  slope  nearshore,  the  water  level  of  the  entire  shelf 
is  below  its  initial  value. 

Figure  3-27  shows  the  juxtaposition  of  sea  surface  height  “snapshots”  every  5  time  steps. 
Notice  that  the  perspective  has  been  chosen  so  that  the  steady-state  surface  profile  can  be  more 
easily  seen.  The  cusp-like  nature  of  the  signal  is  clearly  visible.  It  is  also  easy  to  see  the  rapid 
drop  in  sea  surface  height  and  subsequent  oscillations  after  the  initial  set-up  .  Higher  frequency 
small  scale  oscillations,  which  were  not  apparent  in  figures  3-25a  and  3-26,  are  seen  on  the  face 
of  the  larger  signal.  Similar  motions  appeared  in  the  flat  bottom  cases. 

The  fact  that  the  sea  surface  has  a  minimum  in  the  interior  of  the  domain  is  apparent  as 
steady-state  is  approached.  Clearly  the  cross-shore  momentum  balance  varies  in  y  since  the 
cross-shore  pressure  gradient  differs  so  greatly  across  the  shelf  (even  changing  sign).  Remember 
that  this  did  not  happen  in  the  flat  bottom  case  where  the  steady-state  cross-shore  pressure 
gradient  was  constant  everywhere.  Further  discussion  of  the  pressure  gradient  is  provided  in 
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Figure  3-25:  Time  series  of  a)  sea  surface  height,  b)  depth-averaged  cross-shore  velocity  and  c) 
depth-averaged  alongshore  velocity  at  two  points  on  the  shelf.  The  offshore  point  (solid  line)  is 
located  in  60m  of  water  and  the  onshore  point  (dashed  line)  is  located  in  20m  of  water. 
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Figure  3-26:  Contour  plot  of  sea  surface  height  on  the  y-t  plane.  Tic  marks  are  at  intervals  of 
5A t  and  Ax  on  the  horizontal  and  vertical  axes,  respectively.  The  time  axis  ranges  from  0  to 
40  hrs. 
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Figure  3-27:  3-D  plot  of  sea  surface  height  in  y-t  space.  Tic  marks  are  at  intervals  of  5A t  and 
Ax  on  the  t  and  y  axes,  respectively.  The  time  axis  ranges  from  0  to  40  hrs. 
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the  momenum  balance  subsection  below. 

The  cross-shore  velocity  is  shown  in  figure  3- 25b.  As  with  the  surface  height,  the  oscillation 
magnitudes  are  comparable  to  the  flat  bottom  shallow  case  while  the  decay  time  is  similar  to  the 
flat  bottom  intermediate  case.  Initially,  the  flow  accelerates  faster  at  the  onshore  point.  This  is 
understandable  since,  in  the  absense  of  bottom  stress  and  a  pressure  gradient,  the  cross-shore 
momentum  balance  reduces  to  dV/dt  =  Tay/p  yielding  essentially  equal  transports  across  the 
shelf.  Equal  transports  mean  larger  depth-averaged  velocities  in  shallow  water.  Each  series 
peaks  due  to  the  establishment  of  an  adverse  pressure  gradient  upon  the  arrival  of  the  sea 
surface  signal.  The  flow  is  then  rapidly  accelerated  in  the  opposite  direction  and  soon  turns 
offshore.  The  cross-shore  flow  oscillates  a  few  times  before  reaching  steady-state.  Both  time 
series  show  a  net  offshore  flow  (as  would  be  expected  everywhere  in  light  of  figures  3-26  and 
3-27  having  only  negative  sea  surface  values  at  steady-state). 

Figure  3-25c  shows  time  series  of  alongshore  velocity.  Both  points  spin-up  with  maximum 
speeds  similar  to  the  shallow  flat  bottom  case.  The  difference  between  the  flat  and  linearly- 
sloping  bottom  cases,  however,  is  that  the  steady-state  alongshore  flow  varies  across  the  shelf. 
This  is  entirely  due  to  the  fact  that  zo/h  varies  across  the  shelf  causing  the  DRM  to  produce 
a  variable  frictional  relationship  across  the  shelf.  Cross-shore  variation  cannot  be  produced 
by  a  model  with  spatially-constant  r  (or  Cj)  since  the  steady-state  depth- averaged  alongshore 
momentum  balance  would  not  vary  across  the  shelf  {raxjp  =  rU/h  =  rxalongshore  depth- 
averaged  velocity)). 

3.2.3  Momentum  balances 

The  addition  of  cross-shore  depth  variation  in  this  scenario  implies  that  the  momentum  balance 
at  different  points  across  the  shelf  varies.  Figures  3-28  and  3-29  show  momentum  balances  at 
the  onshore  and  offshore  points,  respectively. 

At  the  onshore  point,  the  alongshore  momentum  balance  (figure  3- 28a)  initially  behaves  like 
a  forced  inertial  oscillation,  but  only  for  the  short  time  that  the  bottom  stress  is  zero.  Once 
the  bottom  stress  is  updated,  it  rapidly  grows  until  essentially  equaling  the  surface  stress.  This 
takes  slightly  less  than  5  hrs  (a  time  comparable  to  the  shallow  flat  bottom  spin-up  time).  The 
approach  to  steady-state  is  subsequently  modulated  by  damped  modified  inertial  oscillations 
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Figure  3-28:  Time  series  of  momentum  balance  components  at  the  onshore  point. 
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Figure  3-29:  Time  series  of  momentum  balance  components  at  the  offshore  point. 


superimposed  on  the  basic  Couette  balance.  The  oscillations  persist  much  longer  than  they 
do  in  the  flat  bottom  shallow  case  indicating  that  waves  are  damped  more  slowly  here.  This 
makes  sense  because  the  waves  spend  some  time  propagating  through  deep  water  where  friction 
is  weaker. 

The  alongshore  momentum  balance  at  the  offshore  point  (figure  3- 29a)  develops  more  slowly. 
It  still  evolves  from  the  initial  forced  inertial  wave  balance  to  a  modulated  Couette  balance, 
then  to  pure  Couette  balance.  However,  the  overall  spin-up  time  is  comparable  to  the  flat 
bottom  intermediate  case.  There  is  a  slight  difference  between  za/h  here  and  in  the  flat  bottom 
intermediate  case,  but  it  does  not  affect  the  spin-up  time  much. 

The  cross-shore  momentum  balances  at  the  inshore  (figure  3-28b)  and  offshore  (figure  3- 29b) 
points  have  many  similarities  with  the  alongshore  balances  described  above.  The  overall  spin-up 
time  at  each  point  is  predictable  from  the  flat  bottom  case  with  the  same  local  value  of  uma/ fh , 
and  the  decay  of  oscillatory  behavior  is  determined  by  the  smallest  ums/ fh  values  in  the  domain 
(i.e.  0(10  hrs)  offshore).  At  steady-state,  the  primary  cross-shore  balance  is  between  the  wind 
stress,  pressure  gradient  and  Coriolis  force  with  bottom  stress  playing  a  secondary  role,  as  it 
did  in  section  3.1.  However,  as  in  the  flat  bottom  scenario,  its  importance  in  modifying  the 
pressure  gradient  field  can  still  have  implications  for  the  surface  set-up. 

3.2.4  Bottom  stress  and  drag  tensor 

In  the  cross-shore  direction,  there  is  a  finite  downwind  steady-state  bottom  stress  at  both 
points  (figures  3-28b  and  3-29b).  Not  surprisingly,  the  cross-shore  bottom  stress  is  bigger  at 
the  onshore  point  than  at  the  offshore  point.  This  implies  the  magnitude  of  the  bottom  stress 
vector  is  slightly  larger  at  the  onshore  point  than  at  the  offshore  point  because  the  alongshore 
bottom  stress  component  is  constant  across  the  shelf.  It  also  implies  that  the  angle  between 
the  bottom  stress  and  depth-averaged  flow  is  slightly  larger  at  the  onshore  point. 

Since  the  steady-state  bottom  stress  magnitude  is  larger  and  the  depth-averaged  speed  is 
smaller  at  the  onshore  point,  the  resistance  coefficient  and  estimated  drag  coefficient  are  larger 
at  the  onshore  point.  Figure  3-30  shows  time  series  of  the  resistance  coefficient  at  the  two 
points.  As  expected,  the  steady-state  values  differ.  The  resistance  coefficient  onshore  is  the 
same  as  it  was  in  the  flat  bottom  shallow  case  since  «*,/ fh  and  zo/h  both  match.  The  resistance 
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coefficient  offshore  is  very  different  from  the  flat  bottom  intermediate  case.  Jenter  and  Madsen 
(1989)  have  shown  that  this  difference  is  due  to  different  roughness  length  to  water  depth  and 
Ekman  depth  to  water  depth  ratios  at  the  two  points. 

There  is  also  a  cross-shelf  difference  in  transient  behavior  of  the  resistance  coefficient.  The 
spin-up  times  and  oscillations  reflect  those  of  the  basic  variables.  The  one  exception  seems  to 
be  an  anomalously  high  value  of  r  at  the  offshore  point  during  early  times.  To  understand  this, 
it  is  useful  to  examine  time  series  of  bottom  stress  magnitude  (figure  3-31)  and  depth-averaged 
current  speed  (figure  3-32)  at  the  offshore  point  from  both  the  DRM  and  the  2|-D  model. 
The  two  figures  show  that  the  large  resistance  coefficient  corresponds  to  a  small  “hump”  in 
the  bottom  stress  concurrent  with  a  relatively  low  velocity.  It  makes  sense  that,  when  the 
depth-averaged  velocity  is  small,  r  (and  even  more  so  Cd)  is  overly  sensitive.  Fortunately,  the 
bottom  stress  is  usually  small  at  these  times  reducing  the  importance  of  the  actual  value  of  r. 

Figure  3-33  shows  the  estimated  drag  coefficient  at  the  two  points  as  a  function  of  time. 
The  same  basic  behavior  as  the  resistance  coefficient  is  evident.  The  anomalous  spike  at  the 
offshore  point  at  early  times  also  appears  in  this  signal,  but  is  even  more  pronounced  because 
Cd  is  roughly  proportional  to  the  inverse  of  depth-averaged  velocity  squared  .  The  steady-state 
value  of  the  drag  coefficient  at  the  onshore  point  is  the  same  as  the  flat  bottom  shallow  case,  but 
the  offshore  value  is  smaller  than  the  flat  bottom  intermediate  case.  The  difference  is  primarily 
due  to  the  difference  in  zo/h  (Jenter  and  Madsen;  1989).  This  is  supported  by  the  fact  that 
there  was  no  variation  in  Cd  among  the  flat  bottom  cases  which  all  had  the  same  roughness 
length  to  water  depth  ratio  (figure  3-18). 

The  angular  behavior  of  the  drag  tensor,  bottom  stress  and  depth-averaged  velocity  is 
depicted  in  figures  3-34  and  3-35.  At  both  points,  the  (3  starts  out  at  45°  and  decreases  linearly 
for  the  short  period  that  there  is  no  bottom  friction  or  pressure  gradient.  After  these  terms 
develop,  it  drops  more  rapidly  until  the  flow  is  approximately  alongshore.  It  then  oscillates 
with  the  cross-shore  velocity  until  reaching  exactly  0°  at  steady-state. 

A  rapid  approach  to  constant  6  can  be  seen  in  both  figures.  Long  before  a  and  (3  reach 
steady-state,  their  difference  does.  This  is  especially  true  in  figure  3-34  which  has  important 
implications  for  using  the  2^-D  model  to  model  time-dependent  forcing  problems.  The  rapid 
angular  response  of  the  drag  tensor  means  that  it  will  not  place  undue  restrictions  on  the  2|-D 
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Figure  3-30:  Time  series  of  resistance  coefficient  at  the  onshore  and  offshore  points 


Bottom  stress  magnitude  —  offshore 


Figure  3-31:  Time  series  of  kinematic  bottom  stress  magnitude,  at  the  offshore  point. 

Solid  line  indicates  DAM  values  and  dashed  line  indicates  DRM  values. 
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Depth-averaged  velocity  magnitude  —  offshore 


Figure  3-32:  Time  series  of  depth-averaged  velocity  magnitude,  \U\fh,  at  the  offshore  point. 
Solid  line  indicates  DAM  values  and  dashed  line  indicates  DRM  values. 
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xlO-3  Estimated  drag  coefficient,  Cd 


Figure  3-33:  Time  series  of  estimated  drag  coefficient  at  the  offshore  (solid  line)  and  onshore 
(dashed  line)  points.  The  spike  in  the  offshore  point  signal  is  truncated  at  .004  although  its 
actual  maximum  is  Cd  «  .023. 
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model’s  ability  to  respond  to  slow  temporal  variations  in  forcing. 

With  variable  depth,  there  is  now  a  difference  among  steady- state  9  values  across  the  shelf. 
Onshore  the  value  is  roughly  —6°  as  it  was  in  the  flat  bottom  shallow  case.  Offshore  it  is  closer 
to  —3°  which  is  approximately  the  flat  bottom  intermediate  value.  The  difference  in  9  across 
the  shelf  was  pointed  out  earlier  in  the  momentum  balance  discussions. 

3.2.5  Differences  betweeen  traditional  and  2|-D  model  results 

Results  from  the  flat  bottom  cases  showed  that  failure  to  include  time-dependence  and  di¬ 
rectional  variability  in  a  drag  law  leads  to  differences  in  both  spin-up  and  steady-state  of  a 
wind-driven  shelf.  In  this  scenario,  the  question  of  spatial  variability  in  the  drag  law  was 
added. 

Numerous  tests  (not  reported  here)  have  confirmed  that  similar  conclusions  can  be  drawn 
regarding  temporal  variability  of  the  drag  law  for  both  the  flat  bottom  and  sloping  bottom 
cases.  Namely,  that  spin-up  is  faster  in  the  2  f-D  model  with  a  variable  drag  tensor  provided 
its  magnitude  remains  less  than  the  constant  resistance  coefficient  used  in  the  traditional  model. 

Figures  3-36-3-41  specifically  address  the  failure  to  include  spatial  variability  in  traditional 
drag  laws.  This  was  not  a  factor  in  the  first  scenario  because  of  the  spatial  uniformity  of  the 
topography  and  the  forcing.  The  figures  correspond  to  wind  angles  of  0°,  30°,  60°,  120°,  150° 
and  180°  respectively.  They  depict  steady-state  surface  profiles  for  different  assumptions  about 
the  magnitude  and  spatial  variation  of  the  resistance  coefficient.  The  four  lines  correspond  to: 

•  the  2|-D  model  result  (solid  with  no  symbols) 

•  the  traditional  model  result  with  r  =  rmjn  from  the  2  |-D  model  (solid  with  x’s) 

•  the  traditional  model  result  with  r  =  rmax  from  the  2  |-D  model  (solid  with  o’s) 

•  a  modified  traditional  model  with  r  =  r(y)  from  the  2|-D  model  (solid  with  *’s) 

Note,  the  last  three  assume  there  is  no  angle  between  the  depth-averaged  velocity  and  the 
bottom  stress  (i.e.  0  =  0). 

It  is  clear  from  the  figures  that  qualitative  behavior  of  the  steady-states  is  not  altered  by 
the  choice  of  bottom  stress  formulation.  However,  quantitative  behavior  is  sometimes  strongly 
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Steady-state  Sea  Surface  Height  —  30  deg.  Wind 


Figure  3-37:  Steady-state  surface  profile  for  a  wind  angle  of  30°.  Curves  are  for  the  2|-D  model 
(  )»  traditional  model  with  r  =  rmjn  (x — x),  traditional  model  with  r  =  rmax  (o — o)  and 
traditional  model  with  r  =  r(y)  (* — *). 


Figure  3-38:  Steady-state  surface  profile  for  a  wind  angle  of  60°.  Curves  are  for  the  2|-D  model 
( — ),  traditional  model  with  r  =  rmjn  (x — x),  traditional  model  with  r  =  rmax  (° — °)  and 
traditional  model  with  r  =  r(y)  (* — *), 


Steady— state  Sea  Surface  Height  —  120  deg.  Wind 


Figure  3-39:  Steady-state  surface  profile  for  a  wind  angle  of  120°.  Curves  are  for  the  2|-D 
model  ( — ),  traditional  model  with  r  =  rTn;n  (x — x),  traditional  model  with  r  =  rmax  (o — o) 
and  traditional  model  with  r  =  r(y)  (* — *). 
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Figure  3-40:  Steady-state  surface  profile  for  a  wind  angle  of  150°.  Curves  are  for  the  2|-D 
model  ( — ),  traditional  model  with  r  =  rm|n  (x— -x),  traditional  model  with  r  =  rmax  (o — o) 
and  traditional  model  with  r  =  r(y)  (* — *). 


Figure  3-41:  Steady-state  surface  profile  for  a  wind  angle  of  180°.  Curves  are  for  the  2|-D 
model  ( — ),  traditional  model  with  r  =  rm;n  (x — x),  traditional  model  with  r  =  rmax  (o — o) 
and  traditional  model  with  r  =  r(y)  (* — *). 


affected.  The  difference  between  curves  in  each  figure  gives  an  indication  of  error  magnitudes 
that  might  be  expected  from  ignoring  0  and  using  a  spatially  constant  resistance  coefficient 
(assuming  it  was  chosen  to  fall  in  the  range  of  values  determined  by  the  2|-D  model  physics!). 

For  this  particular  model  domain,  choosing  a  constant  resistance  coefficient  equal  to  the 
minimum  r  value  from  the  2|-D  model  is  better  than  choosing  r  =  rmax-  Better  agreement  for 
the  small  r  value  is  understandable  since  r  predicted  by  the  2|-D  model  does  not  vary  linearly 
across  the  shelf.  The  onshore  region  in  which  r  «  rmax  is  smaller  than  the  offshore  region  in 
which  r  «  rmjn. 

Since  all  solutions  are  pinned  at  the  offshore  boundary  (i.e.  £  =  0  at  y  =  —L),  the  biggest 
differences  in  surface  height  are  away  from  the  shelf  edge.  Note  that  agreement  of  surface 
height  at  any  point  is  an  integral  effect.  It  implies  correct  reproduction  of  the  average  cross¬ 
shore  pressure  gradient  seaward  of  that  point.  This  explains  why  the  profiles  corresponding  to 
r  =  rmjn  correspond  “better”  with  the  2|-D  model  solutions.  The  curves  which  use  r  =  rmax 
deviate  from  the  2  |-D  results  immediately  inshore  of  the  shelf  edge.  Therefore,  the  integrated 
effect  is  seen  all  the  way  across  the  shelf.  The  fact  that  the  large  resistance  coefficient  is  more 
appropriate  onshore  is  illustrated  by  the  r  =  rmax  curves  being  more  parallel  to  the  2|-D 
model  results  near  the  coast  and  having  maximum  differences  from  the  2  |-D  model  results  in 
the  middle  of  the  shelf. 

It  is  interesting  to  note  that  the  r  =  rmjn  and  r  =  rmax  curves  only  bracket  the  2|-D 
model  curve  for  wind  angles  less  than  90°.  Also,  the  rmax  line  lies  above  the  r  =  rmjn  line  in 
this  range.  For  angles  between  90°  and  180°,  both  constant  r  curves  lie  below  the  2  |-D  model 
curve  and  the  r  =  rmax  line  lies  below  the  r  =  rmjn  line. 

The  reason  that  the  profiles  do  not  bracket  the  2^-D  model  profile  is  that  6  =  0  in  all 
of  the  constant  r  cases,  leaving  an  additional  effect  unaccounted  for.  This  is  apparent  from 
examining  the  fourth  curve  in  each  plot.  This  curve  was  generated  from  a  modified  traditional 
model  where  the  resistance  coefficient  was  allowed  to  vary  in  y  in  exactly  the  same  way  as  the 
steady-state  drag  tensor  magnitude  from  the  2  |-D  model.  For  angles  less  than  90°,  allowing  r 
to  vary  spatially  improves  the  agreement  between  the  2|-D  model  and  the  traditional  model. 
There  is  still  some  disagreement  caused  by  assuming  0  =  0.  In  fact,  the  variable  r  curves 
for  angles  greater  than  90°  show  worse  agreement  with  the  2|-D  model  results  than  do  the 
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r  =  rmjn  curves.  Apparently,  the  two  sources  of  error,  using  constant  r  across  the  shelf  and 
setting  0  =  0,  cancel  each  other  to  some  extent  in  this  range. 

Inversion  of  the  r  =  rmax  and  r  =  rmjn  lines  is  caused  by  Coriolis  accelerations  being 
directed  offshore  for  wind  angles  less  than  90°  and  onshore  for  angles  greater  than  90°.  Since 
larger  r  yields  smaller  alongshore  velocities,  the  cross-shore  Coriolis  force  is  smaller  for  r  = 
rmax-  Therefore,  when  one  wind  direction  is  compared  to  its  compliment,  the  pressure  gradient 
does  not  change  as  much  for  r  =  rmax  as  for  r  =  rmjn,  and  the  r  =  rmjn  curve  jumps  over 
the  r  =  rmax  curve. 

3.3  Linearly- Varying  Depth,  Sinusoidal  Coastline 

This  section  describes  wind-driven  flow  over  a  shelf  whose  depth  varies  linearly  away  from  a 
sinusoidal  coastline.  The  sinusoidal  variation  is  chosen  in  order  to  induce  alongshore  (a:)  as  well 
as  cross-shore  (y)  variations  in  drag  tensor. 

3.3.1  Model  domain  and  wind  stress  field 

The  model  domain  for  this  case  is  essentially  rectangular  (figure  3-1)  with  a  larger  alongshore 
dimension  than  was  used  in  sections  3.1  and  3.2.  The  domain  boundary  is  straight  on  three 
sides,  but  sinusoidal  along  the  coast.  The  domain  extends  62Ax  (one  wavelength  of  the  coastline 
-f  Ax)  in  the  x  direction  and  19.5  Ax  across-shore.  The  amplitude  of  the  coastline  variation  is 
5Ax  where  again  Ax  =  10  km. 

The  bathymetry  varies  as  it  did  in  section  3.2.  At  each  point  along  the  coast,  the  depth 
increases  from  5  m  to  70  m  by  5  m  per  grid  point  in  the  y  direction.  The  amount  of  flat 
area  offshore  (h  =  70  m)  is  determined  by  the  proximity  of  the  end  of  the  linearly-varying 
region  to  the  straight  shelf  edge.  The  shelf  edge  is  positioned  so  that  the  narrowest  part  of  the 
model  domain  has  no  flat  region.  The  bottom  is  covered  by  roughness  elements  corresponding 
to  zq  =  1  cm,  meaning  zo/h  varies  from  2  x  10~3  near  the  coast  to  1.4  X  10-4  offshore  with 
sinusoidal  contours  of  constant  zo/h.  The  wind  stress  is  the  same  as  was  used  in  the  previous 
section  (i.e.  corresponding  to  a  33  m/s  speed  and  45°  direction).  It  is  important  to  note  that 
all  angles  in  this  section  are  measured  relative  to  the  x-axis,  not  a  line  parallel  to  the  local 
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coastline.  Due  to  the  identical  cross-shore  bathymetry  and  wind  speed,  the  model  time  step 
(243  s),  diffusion  time  scales  (315-4375  s),  and  update  frequency  (once  every  5  At)  are  the  same 
as  in  the  linearly- varying  depth  straight  coastline  case. 

One  of  the  most  important  reasons  for  examining  this  scenario  is  that  the  angle  between  the 
wind  and  the  coast  is  not  constant,  but  varies  with  alongshore  position.  Given  the  prescribed 
coastline  configuration,  the  angle  that  the  wind  makes  with  the  coast  varies  from  about  15°  on 
the  leeward  side  of  the  coastline’s  “crest”  to  75°  on  the  windward  side. 

3.3.2  Surface  height  and  depth-averaged  velocity 

Figures  3-42  and  3-43  show  contour  plots  of  sea  surface  height  at  various  times  during  spin- 
up  and  at  steady-state.  The  time  difference  between  the  initial  plots  is  approximately  2.4 
hrs  (35  At).  Therefore,  they  adequately  depict  the  overall  approach  to  steady-state  while 
somewhat  under-resolving  the  seiching.  As  with  the  previous  cases,  the  water  level  within  the 
entire  domain  rises  initially  during  the  propagation  of  the  surface  signal  away  from  the  coast. 
However,  it  does  not  rise  uniformly  along  the  coast  since  the  angle  between  the  wind  and  the 
coast  varies.  The  largest  rise  is  on  the  windward  shore  of  the  “crest”  where  the  wind  angle 
relative  to  the  shore  is  the  largest. 

After  the  initial  rise,  the  sea  surface  height  drops  sharply  as  an  upwind  pressure  gradient 
develops.  The  surface  height  remains  positive  only  within  a  small  area  near  the  windward 
shore,  while  dropping  far  below  its  original  level  in  the  rest  of  the  domain.  Once  this  basic 
configuration  is  reached,  it  persists  to  steady-state.  There  is  still  some  low  amplitude  temporal 
oscillation  and  movement  of  the  bulge  along  the  coast.  This  is  evidenced  by  small  differences 
in  contour  locations  among  the  later  “snapshots”. 

The  steady-state  sea  surface  contains  a  thin  steep  band  of  positive  values  along  the  windward 
side  of  the  shoreline  “crest”.  Here,  the  contours  axe  nearly  perpendicular  to  the  wind  indicating 
that  the  pressure  gradient  is  almost  straight  upwind.  The  leeward  side  of  the  “crest”  shows 
a  large  depression  in  the  sea  surface  extending  well  offshore.  In  shallow  water,  the  pressure 
gradient  is  again  in  the  direction  of  the  wind.  This  is  particularly  true  on  the  left  side  of  the 
depression  where  the  sea  surface  slope  opposes  the  wind.  On  the  other  side  of  the  depression, 
the  pressure  gradient  is  downwind  and  the  contour  lines  are  not  nearly  as  well  aligned.  The 
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Figure  3-42:  Contours  of  surface  elevation  at  a)  2.4,  b)  4.8  and  c)  7.2  hrs. 
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Figure  3-43:  Contours  plots  of  surface  elevation  at  a)  9.6  and  b)  12  hrs  and  c)  steady-state. 
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Figure  3-45:  Depth-averaged  velocity  field  at  4.8  hrs.  Scale:  1  cm  =  60  cm/s. 
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Figure  3-46:  Depth-averaged  velocity  field  at  7.2  hrs.  Scale:  1  cm  =  60  cm/s. 
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Figure  3-48:  Depth-averaged  velocity  field  at  12  hrs.  Scale:  1  cm  =  60  cm/s. 


700 


109 


,  km  from  shelf  break 


pressure  gradient  farther  offshore  from  the  leeward  coast  is  roughly  perpendicular  to  the  wind, 
giving  the  depression  the  bowl-like  structure  visible  in  the  figure.  The  sea  surface  is  constrained 
to  be  0  along  the  shelf  edge  so  the  contours  become  straight  and  parallel  to  the  x-axis  there. 

Figures  3-44-3-49  depict  time  evolution  of  the  depth-averaged  velocity.  Initially,  the  flow 
moves  downwind  as  is  evidenced  by  the  vectors  in  the  lower  lefthand  corner  of  the  first  plot. 
Once  the  surface  signal  propagates  out  from  the  coast,  however,  a  large  change  in  speed  and 
direction  occurs.  The  points  nearer  the  coast  exhibit  an  early  transition  to  fast,  essentially 
shore-parallel  flow.  This  behavior  is  soon  adopted  by  all  the  points  in  the  domain  with  a 
superimposed  oscillation  due  to  the  seiching  motions  of  the  sea  surface.  At  steady-state,  the 
flow  is  largely  shore-parallel  with  larger  velocities  where  the  shelf  is  narrower.  The  somewhat 
strange  looking  velocities  at  the  first  grid  point  in  from  the  coast  are  due  to  both  spatially 
averaging  of  U  and  V  nodes  to  yield  a  velocity  vector  at  the  nearest  £  node  and  to  the  fact 
that  the  finite  difference  grid  does  not  smoothly  approximate  a  true  sinusoid.  These  effects 
disappear  a  short  distance  away  from  the  coast. 

It  is  interesting  to  consider  more  closely  the  direction  of  the  steady-state  flow.  Figure  3-50 
shows  contours  of  (/?).  An  appropriate  question  to  ask,  given  this  information,  is  “  How  nearly 
shore-parallel  is  the  flow?”.  If  it  is  shore- parallel  everywhere  (assuming  smooth,  sinusoidal 
coastline  and  isobaths),  the  contour  lines  should  be  straight  and  parallel  to  the  y-axis.  Also, 
their  values  should  range  between  ±27.5°  for  this  geometry.  Figure  3-50  shows  that  this  picture 
is  accurate  in  a  qualitative  sense,  particularly  along  those  areas  of  the  coastline  which  are  most 
slowly  varying  in  x  (i.e.  at  the  “crest”  and  in  the  “trough”).  When  the  coastline  is  steeper, 
the  picture  of  flow  parallel  to  sinusoidal  isobaths  is  disrupted  by  the  staircase  approximation 
of  the  coastline.  Offshore,  away  from  the  immediate  influence  of  the  coastline,  the  flow  is  well 
described  as  shore-parallel  everywhere. 

3.3.3  Bottom  stress  and  drag  tensor 

Bottom  stress  spins  up  much  like  depth-averaged  velocity.  It  grows  rapidly  near  the  coast 
and  propagates  seaward.  The  stress  becomes  roughly  shore-parallel,  and  then  is  modulated  by 
the  wavelike  adjustment  process.  At  steady-state,  the  stress  is  roughly  shore-parallel  with  the 
largest  stresses  appearing  where  the  shelf  is  narrowest  (figure  3-51).  The  stresses  near  the  coast 
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Figure  3-50:  Contour  plot  of  /?  at  steady-state.  Contour  interval  is  2°. 

reflect  spatial  averaging  and  the  irregularity  of  the  numerical  coastline. 

The  steady-state  resistance  coefficient  (figure  3-52)  varies  from  less  than  0.16  cm/s  to  greater 
than  0.34  cm/s  and  exhibits  a  great  deal  of  spatial  structure.  Clearly  an  assumption  of  spatially 
constant  r  is  inappropriate.  A  bowl-like  structure  is  found  in  the  “trough”  where  the  sea 
surface  was  depressed  and  depth-averaged  velocities  are  at  their  minimum.  Also,  the  maximum 
resistance  coefficients  are  found  offshore  of  the  “crest”  where  the  depth-averaged  velocities  are 
at  their  maximum.  The  apparent  correlation  between  resistance  coefficient  and  depth-averaged 
velocity  suggests  that  it  might  be  interesting  to  estimate  a  steady-state  field  of  drag  coefficients. 
Figure  3-53  shows  contours  of  estimated  Cj.  The  contours  are  nearlyly  sinusoidal  in  the  region 
of  varying  depth  and  roughly  constant  in  the  flat  area  offshore.  This  indicates  that  the  estimated 
steady-state  drag  coefficient  field  is  determined  primarily  by  the  quantities  zo/h  and  u„sf  fh 
and  not  the  angle  between  the  wind  and  the  shoreline  since  the  latter  is  variable  along  isobaths. 
The  results  of  sections  3.1  and  3.2  and  those  of  Jenter  and  Madsen  (1989)  suggest  that  the 
majority  of  variation  in  C &  is  due  to  variation  in  zq/h  rather  than  «*,/ fh. 

Another  important  quantity  to  be  examined  is  the  direction  of  the  bottom  stress.  Figure 
3-54  shows  a  contour  plot  of  a  at  steady-state.  As  with  the  depth-averaged  flow  direction,  the 
contour  lines  are  essentially  parallel  to  the  y-axis  where  expected.  However,  they  are  shifted 
relative  to  the  contours  of  /?.  This  can  easily  be  seen  by  examining  the  0°  contour  on  the 
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-  Figure  3-54:  Contour  plot  of  a  at  steady-state.  Contour  interval  is  2°. 

lefthand  side  of  each  figure.  The  0  =  0°  contour  (figure  3-50)  intersects  the  coast  roughly  in  the 
middle  of  the  “trough”  as  would  be  expected  of  truly  shore-parallel  flow.  The  a  =  0°  contour 
in  figure  3-54,  on  the  other  hand,  intersects  the  coast  to  the  left  of  center  in  the  trough.  The 
a  =  —8°  contour  bisects  the  “trough”.  Therefore,  while  the  approximately  straight  contours  of 
a  parallel  to  the  y-axis  indicate  that  the  direction  of  bottom  stress  is  basically  a  function  of  the 
wind  angle  relative  to  the  coastline,  they  do  not  indicate  that  the  bottom  stress  is  shore-parallel. 

The  difference  between  figures  3-54  and  3-50  is  the  steady-state  field  of  0  shown  in  figure 
3-55.  The  contours  are  surprisingly  smooth  relative  to  those  of  a  and  0.  This  is  even  true  one 
grid  point  from  the  coast,  indicating  that  6  is  a  fairly  robust  quantity.  It  seems  to  be  affected 
only  by  the  “global”  angle  between  the  wind  and  the  coastline  rather  than  the  “local”  wind 
angle  which  varies  dramatically  from  grid  point  to  grid  point. 

Notice  that  the  range  of  0  is  between  —1°  and  —16°  indicating  that  the  traditional  assump¬ 
tion  of  0  =  0  is  not  appropriate.  Also  the  fact  that  it  is  always  negative  means  the  bottom  stress 
is  everywhere  directed  offshore  from  the  depth-averaged  flow,  which  has  important  implications 
for  near  bottom  transport.  The  angular  difference  can  be  seen  by  superimposing  the  vector 
plots  of  steady-state  bottom  stress  and  depth-averaged  flow. 
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Figure  3-55:  Contour  plot  of  6  at  steady-state.  Contour  interval  is  1°. 

3.3.4  Differences  betweeen  traditional  and  2|-D  model  results 

Sections  3.1  and  3.2  have  already  illustrated  many  of  the  differences  between  traditional  and 
2  j-D  models.  Here  a  single  comparison  is  made  and  discussed.  The  traditional  model  is  given  a 
resistance  coefficient,  r  =  0.25  cm/s,  which  is  an  “eyeball”  average  of  the  values  in  figure  3-52. 
Of  course,  0  =  0. 

Differences  in  the  surface  set-up  and  depth-averaged  velocity  fields  will  clearly  be  a  strong 
function  of  the  choice  of  r.  This  is,  in  fact,  common  knowledge  as  many  investigators  view  the 
resistance  coefficient  as  a  tunable  parameter  by  which  they  can  roughly  adjust  model  output 
fields  to  agree  with  measurements.  The  effectiveness  of  varying  r  to  alter  the  steady-state 
surface  height  was  discussed  in  detail  in  section  3.2.  The  same  type  of  variation  in  the  basic 
variables  for  the  sinusoidal  coastline  is  not  addressed  here. 

Directional  differences  in  steady-state  depth-averaged  flow  and  bottom  stress  between  the 
two  models  was  not  addressed  earlier.  The  former  did  pot  exist  in  the  straight  coastline  cases 
and  the  latter  was  somewhat  masked  as  differences  in  Q.  In  the  sinusoidal  coastline  case,  both 
can  differ  between  models.  Figure  3-56  shows  the  steady-state  field  of  (3  from  the  traditional 
model  (a  and  /?  are  the  same  for  the  DAM).  It  looks  very  much  like  the  (i  field  from  the  2|- 
D  model  (figure  3-50).  The  flow  is  roughly  shore-parallel  in  those  regions  of  the  domain  not 
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Figure  3-56:  Contour  plot  of  a  (and,  therefore,  /?)  at  steady-state  from  the  traditional  model. 
Contour  interval  is  2°. 

affected  by  the  jagged  coastline.  In  fact,  the  differences  between  models  are  shown  in  figure 
3-57  to  be  roughly  in  the  range  ±5°.  Variation  in  /?  between  models  is  entirely  due  to  the 
difference  in  bottom  stress  formulations. 

Figure  3-56  also  presents  the  steady-state  bottom  stress  direction  for  the  traditional  model 
since  the  flow  and  stress  are  constrained  to  be  in  the  same  direction.  The  traditional  model 
bottom  stress  is,  therefore,  more  shore-parallel  than  the  2|-D  model  bottom  stress.  Note  that 
the  0°  contour  bisects  the  “trough”  in  figure  3-56. 

The  difference  in  steady-state  a  between  the  two  model  results  is  shown  in  figure  3-58.  The 
values  are  slightly  larger  than  the  differences  in  /?  and  negative  for  the  most  part.  This  means 
that  if  the  bottom  stresses  from  each  model  were  used  to  predict  near-bed  transport  direction, 
the  2|-D  model  would  predict  transport  directed  more  offshore  than  would  the  traditional 
model.  Note  that  the  same  would  be  true  if  the  wind  were  directed  at  135°  as  the  results  from 
section  3.1  presented  in  table  3.1  showed. 
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Figure  3-57:  Contour  plot  of  the  difference  between  steady-state  /?  fields  between  the  2|-D  and 
traditional  models.  Contour  interval  is  1°. 
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Figure  3-58:  Contour  plot  of  the  difference  between  steady-state  a  fields  between  the  2|-D  and 
traditional  models.  Contour  interval  is  2°. 


118 


3.4  Linearly- Varying  Depth,  Straight  Coastline,  Rotating  Wind 

Stress 

The  previous  three  sections  dealt  with  transition  from  rest  to  steady-state  of  an  open  shelf 
forced  by  a  constant  wind  stress.  It  was  demonstrated  that  quantitative  differences  in  both 
transient  and  steady-state  behavior  existed,  but  that  the  steady-state  differences  can  be  con¬ 
trolled  to  some  extent  by  allowing  for  slight  modifications  of  the  traditional  model  such  as 
spatial  (but  not  temporal)  variation  of  the  resistance  coefficient.  In  fact,  introduction  of  a 
correctly  chosen  spatially-varying  drag  tensor  (using  both  r  and  0  equal  to  the  steady-state 
2|-D  model  drag  tensor)  can  lead  to  perfect  agreement  between  the  steady-state  solutions  of 
the  2|-D  and  modified  traditional  model. 

Because  of  this,  it  is  reasonable  to  ask  “Why  is  a  2^-D  model  any  better  than  a  modified 
traditional  model?”  There  are  at  least  three  reasons  why  using  a  temporally-constant  drag 
tensor  field  is  not  as  desirable.  They  are: 

•  Even  though  the  initial  conditions  and  steady-states  might  be  the  same  for  the  two  models, 
the  transient  behavior  will  not  be  the  same,  as  was  shown  in  previous  sections. 

•  It  may  be  very  difficult  to  determine  the  appropriate  drag  tensor  field  a  priori.  Results 
like  those  of  Jenter  and  Madsen  (1989)  can  be  used  for  simple  geometries.  However,  once 
there  is  difficulty  determining  the  direction  of  steady-state  flow  relative  to  the  wind  stress 
(as  in  section  3.3),  there  is  no  reliable  method  for  estimating  the  drag  tensor  field. 

•  Transient  phenomena  such  as  coastal  storms  may  not  have  a  steady-state  response. 

The  first  two  points  were  addressed  earlier.  The  third  is  briefly  discussed  here.  A  simple 
temporally-varying  wind  stress  is  applied  to  an  open  shelf  and  the  basic  variable  and  drag 
tensor  responses  are  examined. 


3.4.1  Model  domain  and  wind  stress  field 


The  simplified  wind  stress  field  is  spatially  invariant  but  rotates  uniformly  in  time.  It  is  de¬ 
scribed  by 
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where  T  is  the  rotational  period  of  the  wind  stress.  The  period  is  chosen  to  be  24  hrs  giving 
the  stress  a  diurnal  periodicity. 

Clearly,  this  is  an  idealized  wind  stress,  but  it  contains  the  essential  features  necessary 
to  evaluate  the  basic  relationship  between  transient  wind  stress  and  the  associated  bottom 
stress.  Actually,  the  level  of  complexity,  is  similar  to  many  studies  of  “storm”-generated  coastal 
circulation  (e.g.  Chao,  1981;  Carton,  1984;  and  Heaps  et  al.,  1988).  There  have  been  more 
realistic  wind  stress  fields  employed  for  studying  wind-driven  coastal  responses,  but  these  have 
usually  been  aimed  at  studying  a  particular  geographic  area’s  response  or  a  specific  type  of 
meteorological  event  (e.g.  Jelesnianski,  1965;  Beardsley  and  Haidvogel,  1981  and  Flather  and 
Khandker,  1987). 

The  magnitude  of  the  wind  stress  corresponds  to  30  m/s  wind  speed  as  used  in  earlier 
scenarios.  Its  direction,  however,  is  linearly  variable  in  time  rather  than  constant.  The  basin 
geometry  is  identical  to  that  in  section  3.2  with  a  straight  coastline  and  a  linearly-sloping 
bottom.  The  bottom  roughness  length  is  again  given  by  Zq  =  1  cm. 

Because  of  the  equivalent  geometries  the  diffusion  and  wave  time  scales  are  the  same  as 
those  in  sections  3.2  and  3.3,  and  that  the  wind  stress  period  is  long  in  comparison.  This  means 
that  the  water  column  has  time  to  respond  to  small  changes  in  wind  direction.  If  the  period 
of  the  wind  stress  were  shorter  than  the  diffusion  time  scale,  the  water  column  could  not  be 
expected  to  respond  appropriately.  As  discussed  in  chapter  2,  the  1-D  model  acts  as  a  low-pass 
filter,  only  communicating  slowly- varying  information  to  the  bottom. 


3.4.2  Sea  surface  height  and  depth-averaged  velocity 

The  response  of  the  shelf  to  periodic  forcing  is  in  many  ways  different  from  the  scenarios 
presented  earlier.  In  particular,  the  flow  field  does  not  reach  a  true  steady-state  (i.e.  one  where 
there  is  no  temporal  change  within  the  domain  after  large  times).  Rather,  it  reaches  a  situation 
where  the  motion  is  periodic.  This  is  clearly  illustrated  in  figure  3-59a-c  where  time  series  of 
surface  height,  alongshore  velocity  and  cross-shore  velocity  are  shown  at  the  same  onshore  and 
offshore  points  as  defined  in  section  3.2. 
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Figure  3-59:  Time  series  of  a)  surface  height,  b)  alongshore  depth-averaged  velocity,  c)  cross¬ 
shore  depth-averaged  velocity  at  the  offshore  (solid  line)  and  onshore  (dashed  line)  points  and 
d)  alongshore  (solid  line)  and  cross-shore  (dashed  line)  kinematic  wind  stress. 
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The  sea  surface  signal  behaves  as  would  be  expected.  The  larger  response  is  onshore  away 
from  the  pinned  shelf  edge.  Also,  the  onshore  point  responds  first  as  the  signal  propagates  out 
from  the  coast.  Consequently,  the  onshore  response  leads  the  offshore  response.  This  can  be 
seen  in  the  alongshore  velocity  signal  too. 

The  initial  response  of  the  sea  surface  at  each  point  shows  a  very  small  rise  followed  by  a 
more  pronounced  drop  because  the  wind  direction  starts  out  at  0°  and  rotates  toward  the  shore 
(figure  3-59d).  There  is  a  small  bulge  pushed  up  against  the  shore  during  the  first  few  time 
steps.  This  creates  the  slight  positive  signal.  Shortly  thereafter,  however,  a  positive  alongshore 
flow  develops  sufficiently  to  move  water  away  from  the  coast  through  Coriolis  acceleration, 
thereby  causing  the  dip.  The  motion  is  nearly  periodic  from  that  time  on. 

There  is  noticable  high  frequency  motion  in  figures  3-59a  and  3-59c.  This  is  due  to  small 
alongshore  pressure  gradients  generated  each  time  the  wind  is  exactly  perpendicular  to  the 
shore.  As  explained  earlier,  the  equations  of  motion  are  degenerate  when  the  wind  is  straight 
onshore.  However,  longer  model  runs  than  those  presented  reveal  that  the  high  frequency 
motions  do  not  grow  from  one  period  to  the  next.  Even  though  the  motions  are  generated 
when  the  wind  is  cross-shore,  they  are  most  strongly  forced  when  it  is  alongshore.  In  figures 
3-59a  and  3-59c,  times  of  maximum  high  frequency  motion  correspond  to  times  when  the  wind 
is  directly  alongshore.  They  are  sufficiently  damped  as  the  wind  rotates  to  render  them  stable 
over  long  times.  The  magnitude  of  the  small  alongshore  pressure  gradients  ranges  from  1  x  10-4 
cm/s2  at  generation  to  1  x  10-3  cm/s2  when  they  are  forced  alongshore.  This  corresponds  to 
alongshore  height  differences  between  grid  points  of  0.1  and  1  cm,  which  are  very  small. 

Figure  3-60  shows  contours  of  surface  height  in  the  y-t  plane.  The  basic  periodic  behavior 
is  immediately  apparent.  There  is  a  slight  indication  of  the  signal  propagating  outward  from 
the  shore  as  is  evidenced  by  the  left  side  of  each  set-up  or  set-down  period  being  cusplike  in 
shape.  Not  surprisingly,  a  cusplike  behavior  was  also  seen  in  section  3.2  (figure  3-26).  There  is 
a  difference,  however,  in  that  the  wind  stress  here  is  periodic  thereby  coupling  the  sea  surface 
response  to  a  particular  frequency.  Therefore,  the  oscillatory  behavior  is  a  combination  of  free 
and  forced  waves. 

The  initial  positive  signal  and  subsequent  drop  propagating  out  from  the  shore  are  seen  on 
the  left  side  of  figure  3-60.  These  are  also  readily  visible  in  figure  3-61  which  shows  the  sea 
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Figure  3-60:  Contour  plot  of  surface  height  on  the  y-t  plane.  Tic  marc  spacing  is  5A t  and 
0.5Ax  on  the  horizontal  and  vertical  axes,  respectively.  Tiine  ranges  from  0  to  50  hrs. 
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surface  behavior  in  three  dimensions.  The  high  frequency  part  of  the  signal  is  seen  as  ripples 
along  the  “peaks”  and  “valleys”  in  figure  3-61.  It  is  somewhat  difficult  to  tell  from  the  figure, 
but  the  high  frequency  signal  is  generated  locally  everywhere  across  the  shelf  at  the  same  time. 
In  other  words,  it  does  not  have  the  signature  of  a  disturbance  created  at  a  particular  cross-shelf 
location  and  propagated  across  the  shelf. 

The  alongshore  velocity  signal  in  figure  3-59b  shows  some  interesting  behavior  when  com¬ 
pared  with  the  results  for  the  linear  depth  constant  wind  stress  case.  Both  scenarios  have  the 
same  geometry  and  wind  stress  magnitude.  The  only  difference  is  that  the  wind  direction  is 
constant  in  section  3.2  and  variable  here.  The  alongshore  velocity  in  the  rotating  wind  case  has 
a  larger  amplitude  onshore  than  offshore.  The  opposite  is  true  at  steady-state  for  the  constant 
wind  stress.  For  short  times,  however,  the  constant  wind  develops  larger  alongshore  flow  at  the 
onshore  point.  This  was  attributed  to  the  more  rapid  response  time  in  shallow  water.  The  same 
feature  appears  in  figure  3-59b.  The  shallow  water  response  time  is  short  enough  so  that  the 
signal  behaves  more  like  a  succession  of  steady-states  while  the  deeper  part  of  the  domain  does 
not  respond  fast  enough.  In  other  words,  the  offshore  point  has  not  approached  steady-state 
by  the  time  the  wind  stress  changes  directions  causing  the  flow  to  seek  a  new  steady-state. 

The  notion  of  slow  response  offshore  is  supported  by  recalling  the  shallow  and  intermediate 
spin-up  times  from  section  3.1.  They  were  5  hrs  and  10  hrs,  respectively.  Clearly,  a  water 
column  that  takes  10  hrs  to  spin-up  is  not  able  to  “keep  up”  with  a  wind  stress  that  varies  with 
a  period  of  24  hrs.  A  water  column  which  takes  5  hrs  to  spin-up  does  a  better  (but  not  perfect) 
job  of  following  the  wind  stress. 

It  is  difficult  to  say  whether  the  amplitude  of  the  alongshore  flow  oscillation  at  the  inshore 
point  in  figure  3-59b  is  equal  to  the  maximum  steady-state  alongshore  flow  due  to  a  constant 
wind  stress,  since  not  all  wind  directions  were  examined  in  section  3.2.  However,  it  is  larger 
than  the  alongshore  velocity  for  a  45°  wind.  The  situation  is  different  at  the  offshore  point.  The 
amplitude  there  is  less  than  the  steady-state  velocity  maximum  for  a  constant  wind,  since  it  is 
less  than  the  constant  45°  wind  result.  This  is  another  indication  that  the  flow  is  not  capable 
of  responding  fast  enough  at  the  offshore  point. 

There  is  another  interesting  phenomenon  in  the  alongshore  velocity  signal.  The  time  series 
is  not  as  sinusoidal  at  the  onshore  point  as  it  is  at  the  offshore  point.  The  main  reason  for  the 
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Figure  3-61:  3-D  plot  of  surface  height  ou  the  y-t  plane.  Tic  marc  spacing  is  5 At  and  0.5 Aar 
on  the  time  and  y  axes,  respectively.  Time  ranges  from  0  to  50  hrs. 


asymmetic  crests  and  troughs  in  the  signal  is  that  sometimes  the  alongshore  wind  stress  and 
cross-shore  velocity  are  in  phase  causing  Tax/p  and  fV  to  have  the  same  sign,  and  sometimes 
they  are  out  of  phase  causing  opposite  signs.  The  former  yields  large  values  of  dU /dt  and, 
therefore,  steeper  variation  in  the  signal.  The  effect  does  not  show  up  offshore  because  the 
phase  relationship  between  the  wind  stress  and  cross-shore  velocity  is  different  there. 

3.4.3  Drag  tensor 

Figure  3-62  depicts  the  time-dependent  behavior  of  the  resistance  coefficient  at  the  onshore 
and  offshore  points.  It  is  clearly  not  constant.  Both  points  start  out  at  0  when  the  domain  is 
at  rest.  The  onshore  point  grows  most  rapidly  at  first  since  the  depth-averaged  flow  develops 
fastest  there.  After  the  initial  spin-up,  the  resistance  coefficient  varies  periodically  with  twice 
the  frequency  of  the  wind  stress.  This  is  because  the  resistance  coefficient  is  positive  definite, 
being  determined  by  the  bottom  stress  and  depth-averaged  velocity  magnitudes.  Each  of  which 
has  two  maxima  during  one  period  of  the  wind  stress. 

Notice  that  the  resistance  coefficient  displays  the  high  frequency  variation  seen  earlier  in 
the  sea  surface  signal  and  cross-shore  velocity.  This  is  understandable  since  the  bottom  stress 
determined  by  the  DRM  is  a  function  of  the  pressure  gradients  in  the  DAM.  Because  the 
high  frequency  motions  indicate  the  presence  of  small  alongshore  pressure  gradients,  the  DRM 
bottom  stress  and  the  resistance  coefficient  are  affected. 

In  addition  to  the  high  frequency  motion,  the  resistance  coefficient  offshore  shows  a  large 
spike  during  each  period,  which  may  be  explained  with  time  series  of  bottom  stress  and  depth- 
averaged  velocity  magnitude  (figure  3-63).  When  the  depth- averaged  speed  drops  to  nearly 
zero,  the  bottom  stress  magnitude  is  small  but  slightly  lagging  the  speed.  Consequently,  the 
resistance  coefficient  is  large  before  dropping  to  almost  zero  as  the  bottom  stress  becomes  small 
and  the  depth-averaged  speed  begins  to  grow.  The  “spike”  phenomenon  was  also  seen  in  figure 
3-30. 

The  fact  that  the  resistance  coefficient  is  not  constant  in  time  leads  to  the  question  “Would 
a  drag  coefficient  be  any  steadier?”  This  was  addressed  to  some  extent  in  section  3.1  where 
it  was  shown  that  the  steady-state  drag  coefficient  was  more  predictable  than  the  steady-state 
resistance  coefficient.  However,  both  quantities  showed  significant  time-dependence.  The  same 
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Figure  3-63:  Time  series  of  a)  depth-averaged  velocity  magnitude  and  b)  kinematic  bottom 
stress  at  the  offshore  point  (solid  line)  and  the  onshore  point  (dashed  line). 
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is  true  for  the  oscillatory  wind  stress  scenario.  Figure  3-64  shows  the  estimated  drag  coefficient 
which  exhibits  as  much  temporal  variability  as  the  resistance  coefficient  3-62. 

As  it  was  for  the  resistance  coefficient,  the  frequency  of  variation  in  Cd  is  twice  that  of  the 
wind  stress.  However,  each  oscillation  of  Cd  is  antisymmetric  about  the  line  Cd  .0035  and  a 
line  half  way  between  each  peak.  Notice  that  the  resistance  coefficient  increases  and  decreases 
nearly  symmetrically  (disregarding  the  peaks  at  the  offshore  point).  This  occurs  because  there 
are  times  when  the  velocity  is  increasing  (decreasing)  more  slowly  than  the  stress  while  at  the 
same  time  the  square  of  the  velocity  is  increasing  (decreasing)  more  rapidly. 

Figure  3-65  shows  the  variation  of  9  with  time,  which  is  clearly  not  constant.  Disregarding 
the  wrap-around  spikes,  there  are  substantial  periods  during  which  9  varies  by  at  least  20°  or 
30°.  It  is  hard  to  imagine  that  a  traditional  drag  law  would  be  appropriate  in  light  of  this 
variation. 

The  temporal  behavior  of  6  is  expected  based  on  previous  results.  Table  3.1  showed  that  for 
the  flat  bottom  shallow  case,  the  steady-state  value  of  9  varied  substantially  with  wind  direction. 
In  fact  the  range  of  values  in  the  table  closely  corresponds  to  the  range  in  the  figure.  Note  also 
that  9  rotated  clockwise  in  table  3.1  as  the  wind  direction  rotated  clockwise,  a  feature  of  figure 
3-65  as  well.  The  unstable  behavior  described  in  section  3.1  for  90°  wind  angles  corresponds  to 
the  wrap  around  here. 

Interestingly,  the  values  in  figure  3-65  are  “90°  out  of  phase”  with  those  in  table  3.1.  Here 
the  wrap  around  occurs  when  the  wind  is  nearly  shore-parallel.  This  happens  because  the 
bottom  stress  and  depth-averaged  velocity  are  nearly  zero  at  those  times,  making  9  extremely 
sensitive  to  small  changes  in  the  bottom  stress  or  depth- averaged  velocity.  For  the  steady  wind 
situation,  depth-averaged  velocities  were  nearly  zero  when  the  wind  was  perpendicular  to  the 
shore. 

Clearly,  the  large  range  of  9  and  its  complicated  phase  structure  make  a  priori  determination 
unlikely.  This  is  merely  one  of  many  strong  indications  presented  in  this  thesis  that  traditional 
drag  laws  must  be  modified  to  include  temporal,  spatial  and  directional  effects  in  order  to 
reproduce  correctly  the  relationship  between  depth-averaged  velocity  and  bottom  stress. 
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3.4.4  Differences  betweeen  traditional  and  2|-D  model  results 

One  of  the  most  interesting  aspects  of  this  scenario  is  the  difference  in  response  times  at  the 
onshore  and  offshore  points.  This  leads  to  differences  in  response  magnitude  between  steady 
and  periodic  wind  cases  and,  in  some  parts  of  the  domain,  to  an  asymmetry  in  the  perioidic 
alongshore  velocity  signal.  The  key  feature  is  the  phase  relationship  between  the  wind  stress 
and  the  flow  field.  As  previous  results  demonstrated,  the  temporal  response  of  any  domain  is  a 
strong  function  of  the  drag  law.  Therefore,  one  would  expect  that  a  traditional  drag  law  would 
significantly  alter  the  phase  relationship  between  the  wind  stress  and  flow  field  in  this  scenario. 

Figure  3-66  shows  time  series  of  alongshore  velocity  at  the  onshore  and  offshore  points 
for  the  2|-D  model  and  for  a  constant  r  model  with  r  =  .25  cm/s,  roughly  halfway  between 
peak  values  for  the  onshore  and  offshore  points  in  the  2|-D  model.  As  expected,  the  phase 
relationship  between  the  wind  stress  and  flow  field  is  significantly  altered.  The  asymmetry  in 
the  signal  at  the  onshore  point  is  gone.  Peaks  for  the  constant  r  case  lag  those  of  the  2|-D 
model.  This  is  the  same  phenomenon  seen  in  section  3.1  when  the  constant  r  case  responded 
more  slowly  to  a  steady  wind  stress. 

Not  surprisingly,  the  magnitude  of  the  alongshore  velocity  is  also  different  between  the 
models.  At  the  onshore  point,  the  traditional  model  over-predicts  the  velocity  because  the 
resistance  coefficient  is  too  small.  At  the  offshore  point,  the  velocity  is  too  small  because  r  is 
too  large.  Clearly,  the  flow  field  is  sensitive  to  the  choice  of  drag  law  in  this  scenario,  as  it  was 
in  previous  sections. 
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Figure  3-66:  Time  series  of  alongshore  velocity  for  the  2|-D  model  (solid  lines)  and  for  the 
constant  r  model  (dashed  lines)  at  the  onshore  point  (large  amplitude  signal)  and  at  the  offshore 
point  (small  amplitude  signal). 
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Chapter  4 


Summary  and  Conclusions 


The  main  objective  of  this  work  was  to  develop  an  improved  means  of  specifying  bottom  stress 
in  wind-driven  depth-averaged  circulation  models.  The  result  is  a  so-called  2^-D  model  which 
uses  a  simplified  1-D  formulation  to  provide  a  drag  tensor  at  each  grid  point  of  a  standard 
depth-averaged  model. 

The  drag  tensor  is  used  to  scale  and  rotate  the  depth-averaged  velocity  to  yield  a  bottom 
stress.  Rotation  eliminates  the  restrictive  assumption  of  traditional  depth- averaged  models  that 
depth-averaged  velocity  and  bottom  stress  are  in  the  same  direction. 

The  results  presented  in  this  thesis  are  intended  to  show  three  things.  They  are 

•  the  role  of  bottom  stress  in  wind-driven  coastal  flows 

•  the  variation  of  the  drag  tensor  in  these  situations 

•  a  comparison  between  the  2|-Dimensional  model  results  and  those  of  a  traditional  depth- 
averaged  model. 

The  major  findings  relevant  to  each  are  summarized  below. 

4.1  The  Role  of  Bottom  Stress 

Perhaps  the  most  important  role  of  bottom  stress  is  in  setting  the  response  time  scale  of  a 
particular  domain  to  a  particular  wind  forcing.  The  two  frictional  parameters  which  determine 
the  response  time  were  shown  to  be 
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•  u*»/ fh  ~  where  is  the  magnitude  of  the  kinematic  wind  stress,  /  is  the  Coriolis 
parameter  and  h  is  the  water  depth. 

•  zo/h  -  where  zQ  is  the  bottom  roughness  length. 

The  typical  range  of  u*s/ fh  for  flows  of  interest  in  this  thesis  is  roughly  5-100,  while  the  range 
of  zo/h  is  10-3-10-9.  The  latter  includes  large  “apparent  roughnesses”  due  to  wave-current 
interaction  and  bedforms. 

There  are  actually  two  time  scales  in  the  wind-driven  response.  The  first  corresponds  to 
a  smoothly-varying  overall  approach  to  steady-state,  and  the  second  to  the  decay  of  higher 
frequency  seiches  that  propagate  back  and  forth  across  the  shelf.  Both  are  functions  of  u»«/ fh 
and  zo/h.  As  was  shown  in  section  3.1,  the  two  time  scales  are  roughly  the  same  in  domains 
of  constant  depth.  The  results  yielded  times  of  5,  10,  and  20  hrs  for  u*s/fh  =  30,  10,  and  5, 
respectively.  Since  the  depth  was  constant  in  all  cases,  zo/h  was  always  the  same.  The  different 
time  scales  were,  therefore,  entirely  due  to  different  values  of  u,s//h.  The  larger  u*s/fh,  the 
shorter  the  time  scale. 

The  two  time  scales  are  different  in  domains  with  variable  depth.  The  overall  adjustment 
time  for  any  particular  point  was  found  to  be  a  function  of  the  local  values  of  u+s/fh  and 
z0/h,  while  the  wave  decay  time  was  determined  by  more  “global”  values.  The  non-local  effect 
makes  sense  because  the  waves  propagate  around  the  domain  spending  time  in  all  water  depths. 
Consequently,  long  decay  times,  comparable  to  those  in  deeper  water,  can  be  experienced  in 
shallow  water.  This  was  demonstrated  in  section  3.2. 

For  time- varying  forcing,  the  ability  of  any  particular  point  in  a  domain  to  respond  to  that 
forcing  is  a  function  of  its  overall  spin-up  time.  If  the  spin-up  time  is  short  relative  to  the  time 
scale  of  the  forcing  then  the  flow  behaves  somewhat  like  a  progression  of  individual  steady- 
states.  This  was  the  case  at  the  onshore  point  in  section  3.4.  If  the  spin-up  time  scale  is  not 
short  enough  then  the  flow  response  is  diminished  because  it  cannot  “keep  up”  with  the  forcing. 
The  offshore  point  in  section  3.4  was  a  good  example  of  this. 

In  addition  to  determining  time  scales,  bottom  stress  also  significantly  affects  the  magnitude 
of  response.  Seiche  and  steady-state  magnitudes  are  related  to  uma/fh  in  much  the  same  way 
that  the  decay  and  overall  spin-up  times  were.  In  domains  of  constant  depth,  the  magnitudes 
are  larger  for  larger  values  of  u»a/ fh.  In  domains  of  variable  depth,  the  steady-state  magnitude 
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is  still  a  function  of  the  local  value  of  u**//h,  but  the  wave  amplitude  is  determined  by  the 
value  of  14*$/  fh  in  the  region  of  its  generation. 

The  results  of  sections  3.1  showed  that  the  bottom  stress  plays  somewhat  different  roles  in 
the  alongshore  and  cross-shore  momentum  balances.  In  the  alongshore  direction,  there  must 
be  a  Couette  balance  at  steady-state  in  which  the  bottom  stress  and  wind  stress  are  equal  and 
opposite.  In  the  cross-shore  direction,  the  primary  steady-state  balance  is  between  the  wind 
stress,  pressure  gradient  and  Coriolis  term.  Bottom  stress  typically  plays  a  secondary  role. 
Nonetheless,  it  was  shown  that  small  cross-shore  bottom  stresses  could  significantly  affect  the 
cross-shore  pressure  gradient  despite  being  a  small  term  in  the  total  balance. 

4.2  Variation  of  the  Drag  Tensor 

The  drag  tensor  was  formulated  so  as  to  have  a  magnitude,  r,  which  is  used  by  the  depth- 
averaged  model  to  scale  the  current  speed  to  the  bottom  stress  magnitude  and  an  angle,  6, 
which  is  used  to  rotate  the  depth-averaged  velocity  to  the  direction  of  the  bottom  stress.  The 
DRM  described  in  section  2.2  was  used  to  estimate  drag  tensors  for  the  DAM.  This  was  done 
by  first  calculating  velocity  profile  and  bottom  stress  estimates.  The  profile  was  then  integrated 
to  yield  a  depth-averaged  velocity  estimate.  Finally,  the  bottom  stress  magnitude  was  divided 
by  the  depth-averaged  velocity  magnitude  to  obtain  r,  while  the  difference  in  their  directions 
yielded  6. 

The  DRM  was  formulated  in  such  a  way  that  it  produced  drag  tensor  estimates  efficiently 
while  retaining  a  great  deal  of  physics.  The  only  simplification  required  to  keep  the  DRM 
quasi-analytical  was  that  the  local  time-derivatives  were  assumed  depth-independent.  This 
limits  temporaily-resolvable  processes  to  those  with  periods  much  greater  than  the  diffusion 
time  scale. 

Comparison  of  bottom  stress  estimates  from  the  DRM  and  DAM  in  section  3.1  showed 
adequate  agreement.  However,  agreement  for  the  large  value  of  u*,///i  was  a  bit  better.  This 
may  have  been  a  reflection  of  the  shorter  update  period  providing  improved  temporal  resolution 
and/or  the  assumption  of  barotropic  acceleration  terms  being  more  appropriate  in  the  shallow 
case.  In  addition,  because  the  update  frequency  was  based  on  an  estimate  of  the  diffusion  time 
scale  and  because  updating  introduced  small  numerical  waves,  separation  of  time  scales  was 
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important  so  that  the  numerical  response  was  not  confused  with  the  physical  response. 

The  temporal  behavior  of  r  and  9  was  reported  extensively  for  three  of  the  four  scenarios 
(sections  3.1,  3.2  and  3.4).  In  all  three,  they  showed  a  great  deal  of  time  dependence.  The  overall 
spin-up  and  wavelike  adjustment  time  scales  for  the  constant  wind  cases  were  comparable  to 
those  of  the  flow  field.  With  a  rotating  wind,  they  exhibited  periodic  behavior  with  twice 
the  frequency  of  the  forcing.  This  is  understandable  since  r  is  a  positive-definite  quantity 
determined  by  only  the  magnitudes  of  the  bottom  stress  and  depth-averaged  flow  which  had 
two  maxima  per  forcing  period. 

The  steady-state  value  of  r  was  found  to  be  a  function  of  «*,/ fh,  zo/h  and  the  angle  between 
the  wind  stress  and  coastline.  The  resistance  coefficient  varied  directly  with  u*s/fh  for  the 
constant  depth  cases  in  section  3.1.  Similarly,  the  largest  values  of  r  were  found  in  shallow 
water  for  the  three  variable  depth  cases.  Although  clearly  related  to  cross-shore  differences  in 
ums/ fh ,  some  of  the  variation  in  sections  3.2-3.4  was  also  due  to  variation  in  zo/h.  The  results 
of  Jenter  and  Madsen  (1989)  support  this.  Variation  of  r  with  wind  direction  was  shown  in 
section  3.3,  but  found  to  be  too  complicated  to  render  it  simply  predictable. 

Steady-state  values  of  9  were,  in  general,  found  to  be  small.  However,  their  effect  on  the 
basic  model  outputs  was  significant.  9  was  shown  to  be  a  function  of  u,a/ fh,  zo/h  and  wind 
angle.  The  most  significant  variation  appeared  to  be  with  wind  angle  (section  3.1).  Rotating 
the  wind  counter-clockwise  corresponded  to  rotating  steady-state  9  clockwise.  9  varied  most 
rapidly  as  the  wind  became  more  perpendicular  to  the  shore.  The  same  sense  of  rotation 
occured  temporally  in  section  3.4.  However,  the  time  of  fastest  rotation  corresponded  to  times 
of  alongshore  winds. 

Estimated  drag  coefficients  were  calculated  by  dividing  the  resistance  coefficient  by  the 
magnitude  of  the  depth-averaged  velocity.  This  was  done  to  examine  the  temporal,  spatial  and 
steady-state  behavior  of  Cd-  The  drag  coefficient  was  found  to  have  as  much  temporal  variability 
as  r,  despite  the  fact  that  r  often  appears  to  be  correlated  with  the  depth- averaged  velocity. 
This  indicates  that  traditional  quadratic  drag  laws  are  not  able  to  reproduce  adequately  the 
transient  relationship  between  bottom  stress  and  depth-averaged  flow. 

Unlike  r,  the  steady-state  value  of  Cd  is  only  a  function  of  zo/h.  Results  from  section  3.1 
showed  that  it  is  not  strongly  dependent  on  «*,/ fh.  It  was  demonstrated  in  section  3.3  that 
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Cd  does  not  depend  significantly  on  wind  angle  either.  This  is  of  some  importance  for  choosing 
a  traditional  quadratic  drag  law  over  a  linear  one,  since  it  may  be  easier  to  estimate  values  of 
zo/h  a  priori  for  a  particular  domain  than  to  estimate  u*s/ fh  and  wind  angles.  This  result 
does  not  exempt  traditional  quadratic  drag  laws  from  problems,  however.  In  addition  to  the 
absence  of  temporal  variability  described  above,  traditional  drag  laws  still  do  not  allow  for  the 
existence  of  an  angle  between  the  depth-averaged  flow  and  bottom  stress. 

4.3  Comparison  of  Drag  Tensor  and  Traditional  Drag  Law 
Results 

Comparisons  betweeen  the  2|-D  model  and  a  traditional  model  showed  both  temporal  and  spa¬ 
tial  differences.  This  is  not  surprising  since  traditional  models  have  temporally-  and  spatially- 
constant  drag  laws  which  restrict  the  bottom  stress  always  to  be  in  the  direction  of  the  depth- 
averaged  flow. 

Temporally,  the  traditional  model,  in  general,  spun  up  more  slowly  from  rest.  This  was  due 
to  the  abnormally  high  friction  at  early  times  which  was  caused  by  the  traditional  model  having 
a  much  larger  initial  resistance  coefficient  than  the  2|-D  model.  Different  spin-up  times  imply 
different  abilities  to  respond  to  temporal  changes  in  forcing,  and,  as  was  shown  in  section  3.4, 
a  domain’s  response  time  greatly  affects  its  basic  variable  fields. 

Spatially,  failure  to  account  for  cross-shore  variation  in  r  was  shown  in  section  3.2  to  restrict 
the  steady-state  alongshore  depth-averaged  flow  to  be  constant  across  the  shelf.  Only  when  r 
was  allowed  to  vary,  could  the  velocity  change.  Similarly,  differences  in  steady-state  flow  field 
magnitude  and  direction  were  demonstrated  in  section  3.3  when  r  was  held  constant. 

Clearly,  improper  specification  of  the  alongshore  flow  field  has  significant  implications  for 
the  cross-shore  momentum  balance  since  the  Coriolis  term  is  so  important  and  the  sea  surface 
slope  is  so  sensitive  to  changes  in  the  alongshore  flow.  This  was  particularly  true  in  this  thesis 
because  the  domains  modelled  were  all  clamped  at  the  offshore  boundary.  Therefore,  the  sea 
surface  response  at  any  point  on  the  shelf  was  an  integral  function  of  all  of  the  surface  slopes 
seaward  of  that  point.  Results  in  section  3.2  showed  that  the  steady-state  surface  profile  was 
strongly  affected  by  the  assumed  cross-shore  variation  in  r. 
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In  addition  to  neglecting  temporal  and  spatial  variations  in  r,  traditional  models  assume 
9  =  0.  This  has  significant  implications  because  small  values  of  6  can  produce  large  changes 
in  model  output,  as  was  demonstrated  in  section  3.1  when  2|-D  and  traditional  models  with 
identical  r  values  showed  20%  differences  in  surface  height  when  6  was  just  —6°.  Likewise,  in 
section  3.2,  the  two  models  showed  significant  differences  in  steady-state  cross-shore  surface 
profiles  even  when  r  was  allowed  to  vary  exactly  the  same  way  for  both  models.  The  reason 
for  this  is  that,  in  these  straight  coastline  cases,  no  predicted  cross-shore  bottom  stress  can  be 
reproduced  by  a  traditional  model  because  V  =  0  at  steady-state. 

4.4  Conclusion  and  Recommendations  for  Future  Work 

The  major  conclusion  of  this  thesis  is  that  the  relationship  between  the  depth-averaged  flow 
and  bottom  stress  is  complicated  and  that  traditional  drag  laws  are  unable  to  reproduce  many 
important  phenomena.  Lack  of  temporal  and  spatial  variability,  as  well  as  restriction  of  the 
flow  and  bottom  stress  to  be  collinear  is  the  primary  problem.  Therefore,  introduction  of  a 
drag  tensor,  calculated  as  it  was  in  this  thesis,  appears  to  be  a  worthwhile  modification  to 
depth-averaged  circulation  models. 

The  results  of  this  thesis  point  toward  areas  for  future  consideration.  Primary  among  these 
is  verification  of  the  2|-D  model  using  actual  measurements  of  wind-driven  coastal  flow  fields. 
The  functional  dependencies  of  the  drag  tensor  indicate  that  proper  specification  of  the  fields  of 
u*s/ fh,  zo/h  and  wind  direction  will  be  critical  for  good  model  agreement.  Consequently,  this 
implies  a  need  for  good  meteorological  forcing  and  bottom  roughness  information.  The  latter 
will  require  development  of  a  formulation  for  Zq  which  allows  it  to  be  determined  by,  among 
other  variables,  the  bottom  type,  bottom  configuration  and  wave  climate. 
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